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ABSTRACT 

In regions of very high dark matter density such as the Galactic centre, the capture and 
annihilation of WIMP dark matter by stars has the potential to significantly alter their evolu- 
tion. We describe the dark stellar evolution code DarkStars, and present a series of detailed 
grids of WIMP-influenced stellar models for main sequence stars. We describe the changes in 
stellar structure and main sequence evolution which occur as a function of the rate of energy 
injection by WIMPs, for masses of 0.3-2.0 M Q and metallicities Z = 0.0003-0.02. We show 
what rates of energy injection can be obtained using realistic orbital parameters for stars at the 
Galactic centre, including detailed consideration of the velocity and density profiles of dark 
matter. Capture and annihilation rates are strongly boosted when stars follow elliptical rather 
than circular orbits. If there is a spike of dark matter induced by the supermassive black hole 
at the Galactic centre, single solar-mass stars following orbits with periods as long as 50 years 
and eccentricities as low as 0.9 could be significantly affected. Binary systems with similar 
periods about the Galactic centre could be affected on even less eccentric orbits. The most 
striking observational effect of this scenario would be the existence of a binary consisting of a 
low-mass protostar and a higher-mass evolved star. The observation of low -mass stars and/or 
binaries on such orbits would either provide a detection of WIMP dark matter, or place strin- 
gent limits on the combination of the WIMP mass, spin-dependent nuclear-scattering cross- 
section, halo density and velocity distribution near the Galactic centre. In some cases, the 
derived limits on the WIMP mass and spin-dependent nuclear-scattering cross-section would 
be of comparable sensitivity to current direct-detection experiments. 

Key words: dark matter, stars: evolution, stars: fundamental parameters, stars: interiors, 
Galaxy: centre, elementary particles 



1 INTRODUCTION 

Observations continue to support the existence of non-baryonic 
dark matter (DM; |Bergstrom||2000| |Bertone et"aT1|2005[ |Clowe| 
|et al.|2 006; Komatsu et al. 20091 with a cosmological abundance 
around five times that of baryonic matter but of unknown composi- 
tion. Weakly-interacting massive particles (WIMPs) are a popular 
and convenient class of dark matter candidates because their weak- 
scale masses and couplings naturally give rise to an appropriate 
thermal relic abundance in the early universe. 

Typical WIMPs, such as the lightest neutralino in supersym- 
metry (Jungman et al. 1996), posses non-zero nuclear-scattering 
and self-annihilation cross-sections. The nuclear-scattering cross- 
section makes it possible for WIMPs to collide elastically with nu- 
clei in massive bodies such as stars, obtain velocities less than the 
local escape velocity and become gravitationally bound ( |Press &| 



|Spergel|1985]|Griest & Seckel|1987||Gould|1987a|b| >. This popula- 
tion of WIMPs will continue to scatter off nuclei in the star, sinking 
down to the core and eventually annihilating with other captured 
WIMPs. 

If enough WIMPs are captured, the structure of the host body 
may be altered by the energy produced in WIMP annihilations, 
or by energy transport caused by the WIMP-nucleus scattering 
events themselves. This was first realised over 30 years ago in 
the context of heavy neutrinos (Steigman et al. 19781. The po- 
tential for transport effects to modify the structure of stellar cores 
was initially developed by |Spergel & P ress ( 1985 1 and Faul kner "&| 
Gilliland ( 1985 1. Implications of annihilation for stellar evolution 
were first explored by Salad & Silk ( 19891 and Bouquet & Salad 
( 1989b). A series of subsequent studies (Gilliland et al. 1986| |Rerv 
zini|1987| |Spergel & Faulkner|1988| |Faulkner & Swenson 
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Bouquet et al. 1989; Bouquet & Salad 1989a, Deluca et al. 
Salati|1990||Dearborn et al.||1990a|b| |Giraud-Heraud et al 



1988 
1989 



1990 



Christensen-Dalsgaard 1992, Faulkner & Swenson 1993) focused 
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upon possible impacts of energy transport by 'cosmion' WIMPs 
designed to solve the solar neutrino problem. With the advent of 
neutrino oscillations this problem has of course disappeared. Fur- 
thermore, the existence of much more stringent limits upon the 
WIMP-nucleon scattering cross-sections (e.g. Desai et al.||2004| 



|Angle et al.|2008l |Ahmed et al.|2009| |Behnke et al.|2008| > and an 
improved understanding of the distribution of dark matter on galac- 
tic scales (e.g. Bertone & Merritt 2005 , Diemand et al. 2007 1 mean 
that the likelihood of seeing changes induced purely by WIMP en- 
ergy transport seems somewhat diminished. Indeed, later efforts to 
constrain WIMP physics via helioseismology have proven fruitless 
( |Bottino et al.|2002| >. 

Recent times have seen a resurgent interest in the impacts of 
WIMPs upon stars, now focussing almost exclusively upon the in- 
fluence of annihilation. Moskalenko & Wai ( 2007]> and |Bertone| 
|& Fair bairn (20081 showed that it could be possible to see white 
dwarfs heated by WIMP annihilation, at the Galactic centre and in 
globular clusters respectively. Spolyar, Freese & Gondolo ( 2008 1 



and |Natarajan, Tan & O'Shea|j2009^ showed that WIMP annihila- 
tion might be able to partially inhibit the formation of PopIII stars, 
resulting in giant, cool, primordial stars supported entirely by anni- 
hilation energy. In previous letters the current authors presented the 
first numerical simulations of the structure and evolution of WIMP- 
burning main sequence stars, employing and comparing both a sim- 
ple static structure code and a preliminary version of the evolution- 
ary code we present here (Fairbairn, Scott & Edsjo 2008; Scott, 
|Edsjo & Fairbairn|2008] >. We found that WIMP annihilation in stel- 
lar cores diminishes nuclear burning and causes them to re-ascend 
the Hayashi track, in agreement with the analytical estimates of 
ISalati & Silk|(T989l >. 

|Iocco| (2008) and |Freese et al.| {2008 1 performed simplified 
capture calculations on models of 'naturally-formed' PopIII stars, 
showing that even if the stars were to form normally, they might 
later accrete sufficient dark matter to alter their appearance. The 
dark matter densities considered in these studies and in that of 
|Spolyar et alT] ( f2008> were confirmed as reasonable by Freese et al. 
P009) , using a more detailed treatment of the collapse of the pri- 
mordial dark matter-gas halo. Both groups went on to consider 
different stages of the pre-main sequence evolution of WIMP- 
influenced PopIII stars: |Freese et al.| |2008 1 employed polytropic 
models in an attempt to understand the evolution of the stars postu- 
lated by |Spolyar et al.H2008J , and |Iocco et~alT1 < |2008) followed the 
evolution from the tip of the Hayashi track using a full stellar evolu- 
tion code. Both found stalling phases, but of different durations and 
at different stages of the stars' formative evolution. |Yoon, Iocco &| 
|Akiyama| ( [2008l l and |Taoso et"aL] ( (2008l l have now presented simu- 
lations of main sequence PopIII stars assumed to have formed nor- 
mally, but then allowed to evolve with the effects of WIMP cap- 
ture and annihilation. The last three studies show extended main 
sequence lifetimes and stalling on the Hayashi track, in agreement 
with our earlier conclusions at non-zero metallicity and the results 
we present here. 

Those working on PopIII stars have referred to WIMP-burning 
stars as 'dark stars', whilst we and others working at non-zero 
metallicities have typically used the terms 'WIMP burners' or 'dark 
matter burners'. In the interests of cohesiveness and simplicity, we 
will simply adopt the former term. We do acknowledge that the 
term 'dark star' is something of a misnomer, since stars burning 
dark matter are not strictly dark. As we shall see in the following 
pages, except for cases where the ambient dark matter density is 
extremely high, their luminosities are at least reduced relative to 
normal stars. 



In this paper, we present a detailed analysis of the effects of 
dark matter capture, annihilation and energy transport upon the 
structure and evolution of main sequence stars, specifically those 
which might exist at the Galactic centre. In Sect. [2] we give the full 
description of the DarkStars code and its input physics alluded to 
in |Fairbairn et al.H2008f and |Scott et all (2008 | l. Sect.|3]presents the 
properties of main sequence dark stars, based upon a grid of stellar 
models covering a range of masses and metallicities. We take up 
the questions of the distribution of dark matter close to the Galactic 
centre in Sect. [4] and the properties of stellar orbits there in Sect. [5] 
In Sect. [6] we present results from further grids of evolutionary 
models computed with realistic treatments of the environment and 
orbits expected near the Galactic centre. We also discuss existing 
and potential observations in Sect. [6] then give some final remarks 
on the prospect of detecting or constraining the nature of dark mat- 
ter through such observations in Sect. [7] 



2 THEORY AND MODELLING 

2.1 Capture, annihilation and energy injection 

The total population of WIMPs N(t) present in a star is given 
( |Jungman et al.| 1 996 ) by the equation 

dN(t) 



dt 



C(t)-2A(t)~ E(t), 



(2.1) 



where C(t) is the rate at which WIMPs are captured, A(t) is the 
rate at which annihilations occur and E(t) is the evaporation rate. 
The factor of 2 in the annihilation term arises because each anni- 
hilation destroys two Majorana WIMPs. In many cases of interest 
evaporation is negligible, but we will return to this point later. 

Many approximations to the full expression for C(t) derived 
by Gould ( 1987a I have appeared in the literature, with widely vary- 
ing accuracies. Here we attempt to present the full theory in a com- 
pact and usable form. We will also build upon the following in 
Sect. [4] when we consider alternative halo models. For a star cap- 
turing WIMPs from an infinitely distant halo, the capture rate is 



C(t) = 4tt 



R , 



f(u) 



wCl v (w) dudr, 



(2.2) 



where r is the local height in the star, u is the incoming WIMP 
velocity before it is influenced by the star's gravitational field and 
f(u) is the WIMP velocity distribution in the halo. The local es- 
cape velocity at a height r is v(r,t), and w = w(u,r,t) = 
y/u 2 + v(r, t) 2 is the velocity an incoming WIMP obtains by the 
time it reaches a height r. Q,~ (w) is the rate at which a WIMP with 
velocity w scatters to a velocity less than v, and thereby becomes 
captured. This formula does not apply to capture from an already- 
bound population of WIMPs, such as occurs in an adiabatically- 
contracting DM-gas cloud. 

For a scattering nucleus of mass m nuc and a WIMP mass m x , 
kinematics dictate that the only collisions able to scatter a WIMP to 
velocities less than v are those where the fraction A of the WIMP 
energy lost in the collision obeys 



^< A< JL 

9 ^ <i n i 

with 



/' - 



H±l 



(2.3) 



(2.4) 
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This is clearly only possible for values of u that obey 



2 ' 

which is equivalent to 

2 , M« 2 
u < — 



(2.5) 



(2.6) 



The partial capture rate is then given by 



Mi I V 4 



F,(A)| 2 dA 



(2.7) 



Here i denotes the ith nuclear species, m is its local number den- 
sity in the star and Fi(A) is the ith nuclear form factor. 9 is the 
Heaviside step function. The total cross-section Oi for scattering of 
WIMPs on the ith nucleus can be approximated as ( Jungman et al. 
[19961 IGondolo et al.|2004b 



I3 2 [a SI A 2 + a SD 1] \(S p ,i) + {Sn^f 



(2.8) 



(2.9) 



where 

(m x + m p ) 
m p (m x + m nuc ) 

is the ratio of the reduced masses of the WIMP-nucleus and WIMP- 
proton systems. Here ogi and <jsd are the hydrogen-normalised 
spin-independent and spin-dependent nuclear-scattering cross- 
sections respectively, Ai is the atomic number of the nucleus, Ji 
is its spin and (S p ,s) and {S n ,ij are the expectation values of the 
spins of its proton and neutron systems, respectively. 
Assuming an exponential form factor 



|F(A)| 2 = exp( 



m x w A 
2E 



(2.10) 



for heavy elements and a delta function for hydrogen, the integral 
in Eq. |2.7| can be performed analytically. Here Eo is the coherence 
energy arising from the characteristic nuclear radius. When m nuc 
is expressed in GeV/c 2 , it can be approximated as 

5.8407 x 10~ 2 

E w : — a GeV. (2.11) 



m n u C (0.91mi£c + 0.3) 2 



Making the further assumption that WIMP velocities in the halo 
follow an isothermal distribution with dispersion v, the velocity 
distribution in the rest frame of the dark matter halo is 



*<•>-£(§) 



3/2 n u 2 

^3 eX P 



3u 2 

2^ 



(2.12) 



/*(«) = fo(u) exp ( 



3v: 
2tf 



with p x the ambient WIMP density. In the frame of a star moving 
with velocity v* through the halo, this becomes 

■ 2 ^ smhjSuv^/v 2 ) 
3uv+/v 2 

Using Eqs. |2.7||2.10| a nd|2.13| it becomes possible to perform the 
velocity integral in Eq. |2.2| analytically. One converts the step func- 
tion in Eq. |2.7| to a finite upper limit of integration w m ax,i(r, t) — 

v(r, t)\fu/ '/Lt-,i and obtains 



C(t) = 4tt 



/ 7-2 [Wi(<W,<(r,t)) " Wi(0)] dr. (2.14) 
Jo 



Here 

Wi(u)= f ^-wQ-^du 
J u 

= W (r, t )^/3 ;| rTK;) 

m x vVi,p,i V 2 1 



(2.15) 



-BGv 2 /(B+G) 



(B + H)-iy(H) 



x e B + 



Sjj, [B(v 2 +v(r,t) 2 )+Hv(r,t) 2 ] 

3 



_ m X rr — ri P d 

2iio juj. 2v 2 

T(X) = T~(X) - T+(X), 

I V-B + X 
for heavier elements, and 

Wh(u) = ^nnu(r,t)p x ^| v(r,tfE + M-,h 



(2.16) 
(2.17) 
(2.18) 



m x vv* V 27r I 2\fB 4/ihB 3 / 2 
x (2a/B[(v* -u)e- B(u+v * )2 + (v* +u) (2.19) 



-B(u+v t ,) 2 +4Bv t ul 



l + 2Bvi\E 



E = V^l erf TVS 



■erf[YB(u + «* 



(2.20) 



for the special case of hydrogen. 

The annihilation rate A(t) is simply the integral of the local 
annihilation rate per unit volume a(r, t) 



A(t) = 4tt f r 2 a(r,t) dr, 
Jo 

which is given by 

1 2 

a (M) = ^^on^r^) . 

Here n x (r, t) is the local WIMP number density in the star and 
(o"af}o is the non-relativistic limit of the velocity-averaged anni- 
hilation cross-section. The energy injected by WIMP annihilations 
e a nn per unit mass of nuclear matter is 

<w(r,t) = V '/ * H oss (r,t , (2.23) 

p*{r,t) 

where /0*(r, t) is the local stellar density and u\ OSB (r,t) accounts 
for the fraction of the WIMPs' rest-mass energy which escapes in 
the form of neutrinos. We assume that the WIMPs annihilate only 
into standard model particles which, apart from the neutrinos, very 
quickly deposit their energy in the surrounding gas. Together with 
the energy injection/removal rate e t rans due to conductive transport 
by WIMPs, this gives the total local WIMP energy term 



(2.21) 



(2.22) 



fWIMp(r, t) = £ a „ n (r, t) + Etrans (r, t). 



(2.24) 



This acts as an additional source term in the standard stellar lumi- 
nosity equation at every height in the star. 

2.2 Conductive energy transport and distribution 

The simplest way to describe the density of WIMPs in a star is to 
assume that they have thermalised with the stellar matter. We as- 
sume that thermalisation occurs instantaneously, as it would be too 
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cumbersome to allow for a non-thermalised evolution of the WIMP 
distribution at the same time as evolving a star (no formalism be- 
yond explicit Monte Carlo simulations exists at this stage for doing 
such a thing). We expect this approximation to break down when 
capture rates change rapidly, such as during the short-term evolu- 
tion of stars on elliptical orbits at the Galactic centre (Sec. |6.2[ l. Our 
primary interest is in the longer-term behaviour of such stars, which 
should not be strongly effected by the thermalisation process. 

Thermalisation can be local, such that the WIMP energies re- 
flect the local temperature at every height in the star, or global, 
such that they reflect only a single overall characteristic tempera- 
ture Tw- In the first case, WIMPs are in local thermodynamic equi- 
librium (LTE) with the stellar matter, whilst in the second case they 
are isofhermally distributed. The isothermal assumption gives the 
Boltzmann distribution for particles in a gravitational well 



n x , lBO (r,t) = N(t) 



fJj 



E-B., /kT 



-m x 0(r,t)/kT w (t) 



(2.25) 



J 4 7rr '2 e -mx'K r '> t )/ kT w(t) dr' ' 



where (j>(r, t) is the local value of the gravitational potential, gj 
is the statistical weight of the jth energy level (1 in this case) 
and Ej — m x (j>(r, i) is its gravitational potential energy. Because 
WIMPs cluster so strongly in the centre of a star, previous analyses 
typically assumed that the area in which they reside can be approx- 
imated by a sphere of uniform density p c (t), the central density. 
This sets Tw(i) to the central temperature T c (i), and gives 



f G Pc (t)r 2 - 



<j>(r,t) 
This takes Eq. |2.25| to 



-3/2 r 3 ' 



with 



r x (t) 



3kT c (t) 

2TYGp c (t)l7 



1/2 



(2.26) 



(2.27) 



(2.28) 



In practice, there is no longer any real reason to prefer Eq. |2.27| to 
Eq. |2.25| The gravitational potential must be tabulated anyway at 
all heights in the star to obtain v(r, t) for the capture calculation, 
and the integral in Eq. |2.25| can be evaluated quickly and easily 
using modern computers. Most importantly, a reasonable estimate 
for Tw (i) can be directly calculated from the structure of the star 
(as we describe below). 

The usefulness of the uniform sphere approximation lies in the 
length scale it defines, r x . This gives the approximate scale height 
of the WIMP distribution in the star, which can be compared with 
the WIMP mean free path 



I(r,t) 



k(r,t) 1 = a ini (r,t) (2.29) 



at the centre of the star to give the Knudsen number of the system, 
1(0, t) 



K(t) = 



»*(*) 



(2.30) 



The Knudsen number indicates whether the WIMPs travel a dis- 
tance less than the scale size on average and transport energy lo- 
cally (K < 1), or typically travel out beyond r x before depositing 
energy non-locally (K > 1). 



In the extreme limit K — > 0, the energy transport is com- 
pletely local. In this case WIMPs scatter about so often that they 
are in LTE with the nuclei, and the energy transport is exactly 
the case of LTE conductive transport by a gas of massive parti- 
cles. The exact form of the energy injection/removal rate at a given 
stellar radius (i.e. the contribution to the local luminosity), and 
the corresponding density structure for LTf-distributed (rather than 
isothermally-distributed) WIMPs has been calculated by Gould & 
|Raffelt| ( [T9^0b) . These are 



Jtrans.LTE 



(r,t) = 4Tvr 2 n(r, t)n x ,LTE(V, t)l(r, t) 



x rk^|^ k d ^M (2.31) 
L m x J dr 



and 



n x ,LTE 



(r 5i ) = n x , LTE (0,*)p^] 3/J 

[■ 



x exp 



(*) 



ka(r', t) +m x — 



d<j>(r',t) 



■ dr 



'] , (2.32) 



kT*(r',t) 

where normalisation to J B * 47rr 2 n Xi LTE(^, i) dr — N(t) dictates 
the value of n x ,LTE(0, i). Noting that in general 



L(r,t) = 4tt 



2 'p(r' \t)e(r ,t) dr' ', 



(2.33) 



(2.34) 



we see that 

, ,n 1 dZ/trans,LTE(V, t) 

e trans ,LTE(r,t) = WpM Ar 

Note that our sign convention differs from that of Gould & Raffelt 
( 1990b); Eqs. |2.3l1 and |2.34] both refer to the energy injection rates, 
not the rate at which WIMPs remove energy (cf. the sign conven- 
tions in Eq. |2.24} , 

The factors a and k are the dimensionless thermal diffusiv- 
ity and conductivity, respectively. These vary throughout the star 
according to the relative abundances of the different atomic nu- 
clei (and the distribution of WIMP-nucleus relative velocities, if 
the scattering cross-section is velocity-dependent due to a vector 
coupling to quarks; this is not the case for the neutralino). They 
are obtained through numerical solution of the Boltzmann collision 
equation for any given gas mixture. Gould & Raffelt ( 1990b I found 
and tabulated the values of a and « for gases consisting of WIMPs 
and one other nucleus, varying the WIMP-to-nucleus mass ratio p 
from to 100. Whilst the rigorous thing to do for each physical 
mixture would be to re-solve the Boltzmann equation with a com- 
posite collisional operator given as a linear combination of the op- 
erators applicable to the single-nucleus case, a good approximation 
is to simply take a weighted mean of the tabulated a and k values 
themselves. That is, 



aiUi(r, t) 



Oii(pi 



and 



K(r,t) = {l(r 1 t)^2[Ki(p)l i (r,t)] '} . 



(2.35) 



(2.36) 



For p < 100, a and k can be found by interpolation in the tables 
of |Gould & Raff elt (1990b). For larger values of p, the authors 
found the limiting behaviour a — > 2.5 and k — > ^^/2np. To 
get a smooth curve for both a and k, we set the final point of the 
interpolation table to the limiting values at p = 120, so that for 
p > 120 no interpolation is required and the analytical limits are 
used. 
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Notice that Z/trans,LTE(-R*, t) is always zero, as 
dT dr"'^ C^*) ~ ; tms re fl ec ts the fact that conduction by 
WIMPs never constitutes a net energy source nor sink in the star, 
and any energy the WIMPs remove from a hotter region is always 
returned in full to a cooler region. As shown in Fig.[T] etrans.LTE 
is negative in the inner part of the star, increases with radius until 
it becomes positive, peaks, then drops again to asymptotically 
approach zero as r — > 7?* . 

The statement that there is no net energy outflow due to 
WIMP-nucleon scatterings is equivalent to there being no evapo- 
ration of WIMPs from the star. That is, a captured WIMP never 
upscatters sufficiently energetically from a nucleus to become un- 
bound and exit the system. That this is the case in the LTE regime 
is immediately apparent: WIMPs do not travel far from the centre 
before re-scattering and sinking back to the core. In the non-local 
regime, things are not so clear; when the WIMPs' mean free paths 
are much longer than the system scale height, evaporation could in 
principle be significant. In this case, the evaporation rate must be 
equivalent to the net outward energy flux due to WIMP conductive 
energy transport in the isothermal (non-local) picture. That is, 



rans, iso 



(R*,t) = m x c 2 E(t) = L cvap (-R*,£) 



(2.37) 



£trans,iso(V, t) = 4n / r 2 p(r,t)e t 



>(r,t,Tw)dr (2.38) 



j(r, t) — 4nm x 



r 2 R(r, t, T w ) dr, 



(2.39) 



where R(r, t, Tw) is the local WIMP evacuation rate per unit vol- 
ume from a shell at height r, discussed in detail by Gould ( 1987b I. 
The theoretical rate of energy conduction in the isothermal regime, 



£trans,Sp(V, t, Tw) = 



8,Mk 3 / 2 



(r,t)[T*(r,t)-T w (t)] 



^ m x m mc ,, fT+{r,t) Tw(t)y/" 

> <T i n l (r,i) 7 — — i H — 

(m x + m nuc ,i)^ V m nuCil m x / 



(2.40) 



was developed by Spergel & Press ( 1985 1. The equality in Eq. |2.37| 
uniquely determines Tw, so given a a value of E(t) or a particular 
choice of function R(r, t, Tw), one simply implements an appro- 
priate root- finding algorithm to get Tw- Knowing Tw, one has all 
the information required to compute the isothermal density distri- 
bution (Eq. |2.25] l. Whilst the boundary condition itself (Eq. |2.37) 
arises from the consideration of evaporation, and is important to in- 
clude in order to obtain Tw, in most cases of interest evaporation 
is negligible so E(t) = 0. 

Treating conductive energy transport by WIMPs when if > 1, 
where the LTE conduction approximation breaks down entirely, 
is in general rather difficult. In a second paper, Gou ld & Raffelt| 
(1990a) showed via explicit Monte Carlo solutions to the Boltz- 
mann equation that there is no good way of analytically determin- 
ing the WIMP energy transport for large K. One would naively 
expect that as K — > oo, the WIMPs would behave essentially 
isothermally, following the isothermal density structure (Eq. |2.25fr 
and transporting energy according to Eq. |2.40| Whilst they do get 
rather close to the isothermal density structure, the energy trans- 
port via Eq. |2.40| cannot be reconciled with the true Monte Carlo- 
derived ew in any systematic way. Eq. |2.40| should therefore not 
be used as a description of WIMP conductive energy transport in 
practice, even when K — + oo. Its value is in providing a means of 
treating evaporation which is consistent with the isothermal density 
distribution, such that the characteristic temperature of the distri- 
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Figure 1. A snapshot of the energy deposited by WIMP conductive en- 
ergy transport with height in an example star (1 Mq, Z = 0.01, evolved 
in an isothermal halo with p x = 10 10 GeVcm - 3 , v+ = 220 km s -1 and 
v = 270 km s~ 1 ). The snapshot was taken early in the star's evolution, be- 
fore its structure had time to significantly adjust to the effects of WIMP an- 
nihilation. Conduction by WIMP-nucleus scattering consumes energy from 
the very centre of the core and redeposits it further out. The quantity shown 
here is the final form of etrans (Eq. |2.43| , but it d iff ers fro m e t rans,LTE 
only by some correction factors given in Eqs.|2.4l1and|2.44] 



bution naturally results in the correct evaporation rate (even if that 
happens to be zero). 

In their earlier paper, Gould & Raffelt ( 1990b I showed that as 
K increases, the true conductive luminosity is suppressed relative 
to their analytical prediction (Eq. |2,3l| l. The breakdown occurs in a 
clearly quantifiable way for increasing K, so one way to treat con- 
ductive energy transport by WIMPs in the non-local regime is to 
adopt the local expression, but with a 'semi-empirical' luminosity 
suppression pre-factor in line with the suppression seen numeri- 
cally. The suppression shows a sigmoidal shape on a log K scale in 
|Gould & R affelt s results, 

1 1 



i 



1 _|_ e -(lnJf-hK )/r 



. K 



(2.41) 



The relaxation scale must be about 0.4-0.5 to fit the numerical 
result well; with r = 0.5 this function agrees exactly with the 
suppression function chosen by Bottino et al. ( 2002), so we use 
this value. Ko is the 'crossing point' from the local to non-local 
regimes, where WIMP energy transport is most effective; |Gould &] 
Raffelt|found this to be K « 0.4. 

|Gould & Raffelt (1990b) also showed a similar suppression 
of the LTE WIMP conductive luminosity with radius. This can be 
factored into the final expression for Lt{t, t) as a further multi- 
plicative suppression factor h(r, t), such that 



£trans,ftoal(r > *) = f(K)ff(r,t)L t 



rans. LTE 



(r,t). 



(2.42) 



The final form of the term appearing in Eq. |2.24| and the stellar 
luminosity equation becomes 



^trans — 



4nr 2 p(r, t) dr 



[f(KMr,t)L 



trans, LTE 



(r,t)], (2.43) 



in analogy with Eq. |2.34| The radial suppression function appears 
as a simple cubic polynomial 



f)(r,t) 



<(*) 



1 



(2.44) 
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in lGo uld & Raffelt's numerical results. 

For any given K, the degree to which Ltrans(f) is suppressed 
relative to the LTE prediction can be thought of as a general indi- 
cation of the goodness of the LTE approximation. Since the density 
structure in the extremely non-local regime is at least moderately 
well-described by the isothermal approximation (Eq. |2.25| l, the 
'goodness-of-LTE-assumption' function f(K) can also be used as 
a way of interpolating between the isothermal and LTE (Eq. \2.32\ 
densities to give a semi-realistic density structure for all K: 

n x ,flnai(r,t) = f(if)n x , LTE + [1 - f(K)]n x , iao (2.45) 

Since ^ f(K) 1 and both the LTE and isothermal densities 
are individually normalised, n X! fl na i(r, t) will stay correctly nor- 
malised for all K. Since it is only n Xi LTE which enters into the 
expression for the final conductive energy transport, the density 
structure given by n Xi ft na i is important only for determining the 
radial distribution of injected annihilation energy. 

Because no simple measure of the overall efficiency of con- 
ductive energy transport by WIMPs presently exists in the litera- 
ture, we define such a quantity as 



<£(*) = 



2/3*0", t) 



Ctr 



H*(r,t) Cpthcr 

dr 



dr 



K * 2 p*(r,t) 



(2.46) 



the dimensionless WIMP conductive effectiveness. Here (U* (r, t) is 
the mean particle weight, e ther refers to all other energy terms (nu- 
clear, gravitational and annihilation), and the denominator is simply 
a normalisation factor. 8*(t) is therefore the volume- and number- 
density- weighted, integrated ratio of WIMP-mediated energy trans- 
port to all other energy terms. The weighting by number density 
is appropriate because WIMP energy transport is presumably more 
relevant in areas of higher nuclear density. The absolute value arises 
because etrans is a transport rather than a net source term (i.e. takes 
on both positive and negative values). Roughly speaking, in a star 
where WIMP conductive energy transport is the most important 
local source of luminosity, £ is greater than 1 ; where it is a sub- 
dominant contributor, l£ is less than 1. Whilst the weighting with 
local number density rightly biases <E(t) towards the stellar core, 
this means that it is possible for extremely effective, localised en- 
ergy transport by WIMPs in the very centre of a star to dominate the 
overall stellar energetic effectiveness, without significantly altering 
the overall structure. This is because in such a case, the enhanced 
transport only occurs in the most central parts of the core. 



2.3 The DarkStars code 

DarkStars includes WIMP capture based on Eq. | 2.2| generalised 
from the solar capture routines of DarkSUSY (Gondolo et al. 
|2004| >. Exponential form-factor suppression (Eq. |2. 10 1 is assumed 
for scattering off nuclei heavier than hydrogen. In the basic ver- 
sion, the integral over incoming WIMP velocities can either be 
performed analytically assuming an isothermal velocity distribu- 
tion (using Eqs. |2. 15| and |2.19 1, or numerically over any arbitrary 

■ we describe two additional velocity 



velocity distribution. In Sect.H 1 
distributions which we have implemented in the code, one of which 
includes another analytical option for the velocity integral. 

The capture routines are coupled to the EZ version (jPax-| 



ton 



2004|> of the STARS stellar evolution code ( |Eggleton||1971| 



1972| Pols et al.||l995[ >, which uses relaxation to solve the hy 



radial mesh at each timestep. We implement annihilation accord 
ing to Eqs .|2.21f}2.23| W e determine the local WIMP density with 
Eqs. |2,25)|2.32| and |2.45| obtaining T\y as the solution to Eq "~ " ~ 
We include conductive energy transport via Eqs. |2.3l] and Eq 



2.37 



2.43 



The WIMP population is advanced at each timestep by solv- 
ing Eq. |2.1| We assume that since the stellar structure changes very 
slowly under the influence of the WIMPs in comparison to the evo- 
lution of N(t), the evolution of the population between timesteps 
can be well described by the solution to Eq. |2.1| in the special case 
where C(t) and A c (t) = are constant in time: 



N(t + At) = < 



Teq(t) 



C(t)r e<1 (t) 

N(t) 



tanh I 



At 



i 

(t)+t cquiv 



l + N(t)A c (t)At 



VC(t)A c (t)' 



^equiv — t&nh 



. f dN(t) \ 

5 = Blsa {-^r)- 



N(t) 



C(t)r cq (t) 



when C(t) / 0, 

when C(t) = 0, 
(2.47) 
(2.48) 

(2.49) 
(2.50) 



drostatic equations of stellar structure over a 199-point adaptive 



This is well-justified because both C(t) and A c (t) depend only 
upon the stellar structure, not directly upon the absolute WIMP 
population. Here T cq (t) is the emergent time-scale of equilibration 
between capture and annihilation, and i cqu iv is an equivalent earlier 
time from which the approximate solution needs to be evolved for 
the current values of C and A c , Our approximation is an example of 
the general approach to solving stiff differential equations by sepa- 
ration into fast and slow subsystems known as coarse-graining, and 
allows a numerical solution to Eq. |2.1| with timesteps of the order 
of those typically required for stellar evolution. 

This scheme constitutes an explicit solution to Eq. |2.1| where 
each new stellar model is converged with the WIMP population at 
the previous timestep, which is calculated with capture and anni- 
hilation rates computed using the stellar structure of the previous 
model. The models are therefore not completely self-consistent, as 
the WIMP population lags the stellar structure by one timestep. 
Implementing Eq. |2.1| in the internal implicit differencing scheme 
of the STARS code would have required extensive revision of the 
internal solver. As a consistency check, we have implemented a 
'reconvergence mode' similar to that described by Dearborn et al. 
( 1990b), where models are reconverged with the new WIMP pop- 
ulation at every timestep, producing a fully self-consistent solu- 
tion. We have also experimented with rescaling the automatically- 
chosen timesteps to smaller values. Except for some special cases 
which we describe in Sect. [3] the results do not change. 

In general we limit timesteps to allow no more than a certain 
proportional change in the WIMP population per step. Typically we 
demand that the population does not change by more than the cur- 
rent value in one step, but for the more extreme situations in Sect. [6] 
we reduce this by a factor of ten. To aid initial convergence and pre- 
vent this prescription from demanding impossibly small timesteps 

early in the simulation, we begin simulations with populations of 
10 30-35 WIMPs . This 

is many orders of magnitude less than the 
population required to have an effect upon the stellar structure. 

We calculate capture by the 22 most relevant nuclei: X H, 3 He, 
4 He, 12 C, 13 C, 14 N, le O, ls O, 20 Ne, 24 Mg, 23 Na, 27 A1, 28 Si, 32 S, 
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Figure 2. Initial capture rates (top) and annihilation-to-nuclear luminosity 
ratios (bottom) achieved by stars evolved in differing dark matter densities. 
Nuclear luminosities L nuc (0) are initial values, whilst annihilation lumi- 
nosities Lw.max are the maximum values achieved during a star's life- 
time. The dark matter velocity structure is that typically assumed by direct 
detection experiments, where the distribution is isothermal with solar val- 
ues v* = 220 and v = 270 kms . We refer to this canonical example 
as the 'reference solar configuration' (RSC), where it should be understood 
that only the velocity structure, not the density, is reflective of the true so- 
lar situation. These plots provide a simple conversion mechanism between 
capture rates, WIMP luminosities and equivalent RSC dark matter densities. 



3 Fe, 58 Ni, Ni, 



J Pb, 



7 Pb and 207 Pb. The stel- 



lar code follows the abundances of J H, 4 He, 12 C, 14 N, 16 0, 20 Ne 
and 24 Mg, and we assume that the remaining mass is distributed 
amongst the other 15 species according to their abundance ratios 
in the Sun. The data on inter-elemental ratios comes from IAs-1 
plund, Grevesse & Sauval ( 2005 , but with the Ni abundance from 
Scott et aLj 2009 ), and on iso topic ratios from |Heber et al.| ( 2003 
a He/ 4 He), |Scott et al.||2006l 12 C/ 13 C and 16 Q/ 18 0), preliminary 
results from work in progress ( 58 N i/ 60 Ni) and 



Tatsumoto et al 



( 1976| via |Anders & Grevesse|l989| 208 Pb/ 207 Pb/ 2Ub Pb). Atomic 



and nuclear masses are sourced from Wapstra et al. ( 2003]l and |Audi| 
|et al.| {2003} , The present code allows total heavy-element mass 
fractions Z of 0.0001-0.03, which are paired with corresponding 
helium mass fractions of 0.24-0.30. 

Since the annihilation rate goes as n 2 whilst the evaporation 



Figure 3. WIMP-to-nuclear luminosity ratios achieved my a 1 Mq star with 
different annihilation cross-sections (top) and WIMP masses (bottom). The 
dark matter halo configuration is the RSC. Since capture rates do not depend 
upon the annihilation cross-section, the WIMP luminosity is independent of 
the annihilation cross-section in stars where capture and annihilation are in 
equilibrium at the time of maximum annihilation luminosity. Because 
and v are relatively large in the RSC, WIMP luminosities show a strong 
dependence upon the WIMP mass. 



rate goes only as n x , the evaporation rate is much smaller for the 
high ambient WIMP densities and capture rates we are typically 
interested in. Even for the Sun, Gould (1987b I found that evap- 
oration is insignificant unless the WIMP mass happens to be rel- 
atively closely matched with a nucleus found there in significant 
abundance. The heaviest such element is iron, so evaporation can 
be considered negligible in the Sun for m x > 60 GeV, which is the 
case for most WIMPs considered interesting today. (The limit given 
o x > 4 GeV, but this assumed that elements heav- 



by Gould was m x 
ier than helium could be neglected because, at the time, a WIMP 
mass higher than ~ 10 GeV was considered unlikely). We therefore 
simply obtain Tw with E(t) = 0, neglecting evaporation. 

To estimate the factor vi OSE in Eq. |2.23| we carried out ex- 
plicit Monte Carlo simulations of WIMP annihilation in the Sun, 
along the lines of Ble nnow et al.| ( [2"0"0"8"] >. We considered a range of 
masses and annihilation channels, and included full three-flavour 
neutrino oscillations, neutrino interactions, stopping of muons and 
interactions of heavy mesons (e.g. B mesons) in the Sun's core. For 
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essentially all annihilation channels except t + t~, is\ obb for a 100 
GeV WIMP is 5-15% of the rest-mass energy. For heavier WIMPs, 
i^ioss is reduced due to neutrino interactions in the star. For annihi- 
lation to t+t~, v loBB ~ 35-40% for a 100 GeV WIMP (and drops 
for heavier WIMPs). Since the neutralino (our canonical WIMP) 
has very limited annihilation to t + t~, we assume a flat neutrino 
energy loss of v\ OBB = 10% for all annihilations. We also neglect 
the slight dependence upon stellar structure and mass. 

Capture integrals in DarkStars are performed with QUAD- 
PACK C Piessens et al.|1 983 1, whilst simple integrals are done by the 
users' choice of Simpson's rule, Romberg integration or fifth-order 
Runge-Kutta with an adaptive step size. Sufficiently well-behaved 
functions are interpolated using cubic splines. For the others, we 
found the tensional spline routines of Renka (1993) and Testa & 
|Renka| (1999) , after a slight readjustment of the convergence pa- 
rameters, invaluable. 

Some provision has been made in the code for later allowing 
R(r, t, Tw) 7^ 0, alternative form factors and metal-free evolu- 
tion if required. DarkStars is available for public download from 
|http : / /www . f ysik . su . se/~pat/darkstars/| 

Except where explicitly stated otherwise, we perform all sim- 
ulations in this paper with a canonical WIMP mass of m x — 
100 GeV and an annihilation cross-section {a a v)o = 3 x 
10 _26 cm 3 s _1 , which arises from relic density considerations 
assuming WIMPs to be the dominant component of dark mat- 
ter. We use nuclear-scattering cross-sections corresponding to the 
maximally-allowed experimen tal values usi = 10~ 44 cm 2 (Angle 



et al.|2008[|Ahmed et al.|2009|> and <t S d = 10" 38 cm 2 (Desai et al. 



20041 |Behnke et al.|2008) . 



3 BENCHMARK IMPACTS ON MAIN SEQUENCE 
STARS 

To understand the general effects of dark matter accretion and an- 
nihilation upon main-sequence stars, we start by considering a set 
of benchmark stars in a reference halo of WIMPs. In later sec- 
tions, we will expand on these results for more realistic scenar- 
ios. For the benchmark stars, we evolved a grid of models with 
0.3 M < M* < 2M Q , 0.0003 < Z sC 0.02 and 5 < 
< 11. The models were started from the zero- 



lo Sio( GeV cm" 3 i 
age main sequence (ZAMS), and evolved until one of the following 

stopping criteria was met: 

(i) The star left the main sequence (as indicated by the central 
hydrogen mass fraction X c dropping below 10~ 6 ). 

(ii) The star reached a stable equilibrium where all its energy 
was effectively provided by WIMP annihilation (as indicated by 



X c , log 10 p c and log 10 T c changing by less than 10 



-14 M, 



, 10" 



and 10 respectively over four consecutive timesteps). 
(iii) The age of the star exceeded the age of the Universe. 

Except for main-sequence lifetimes, we present results at Z = 0.01 
and just give a brief discussion of the effects of metallicity in the 
text, since most of the properties we discuss did not show any major 
dependency upon metallicity. 

For this grid we used the standard isothermal velocity distri- 
bution, with the default solar values of v* = vq = 220kms~ 
and v — ^/3/2u0 = 270 kms -1 . The resultant capture rates (at 
t = 0) and ratios of annihilation-to-nuclear luminosity are pre- 
sented in Fig. [2] For the sake of comparison, nuclear luminosities 
L nuc are taken at zero age, whilst WIMP annihilation luminosi- 
ties Lw.max are the maximum values the star achieves during its 



evolution. Since stars start with almost no WIMPs, WIMP lumi- 
nosity at zero age is essentially nil, and nuclear luminosity changes 
significantly after zero age as the WIMPs begin to influence the 
stellar structure. As expected, capture rates and WIMP-to-nuclear 
luminosity ratios increase linearly with p x , and lower-mass stars 
capture less but burn a greater ratio of WIMPs to nuclear fuel than 
their higher-mass counterparts. Significantly, WIMP annihilation 
outstrips nuclear burning in a large area of the parameter space. 
Capture rates increase slightly at lower metallicity because of the 
dominance of spin-dependent scattering and capture by hydrogen, 
but are outweighed by the increased nuclear luminosity, causing a 
small decrease in WIMP-to-nuclear burning ratios. 

Because the primary factor governing the impact of WIMPs 
upon stellar evolution is simply N(t), the benchmark results we 
present in this section will hold in general for other combina- 
tions of input particle and halo parameters, subject to an appro- 
priate rescaling. In particular, all scenarios which result in the same 
product of the capture rate and WIMP mass map to roughly the 
same ratio of WIMP-to-nuclear burning, which in turn maps in an 
essentially one-to-one manner to all physical changes in a star's 
structure and evolution. (That is, ignoring 'higher-moment' factors 
like WIMP thermalisation, distribution, conductive energy trans- 
port and capture-annihilation equilibration). In this way, Fig.[2]acts 
as a conversion table between capture rates, WIMP luminosities 
and equivalent dark matter densities in the RSC. (We use the term 
'WIMP luminosity' as shorthand for the ratio ^ w '"q) when it is 
clear what we mean, or when the distinction is irrelevant). 

Lines in Fig. [2] do not all extend to log 10 p x = 11. This is 
because as WIMP luminosity becomes a more significant contribu- 
tor to a star's energy generation, the stellar models become steadily 
more difficult to converge. In many cases, we had to very carefully 
adjust the initial timesteps in order to converge models near the 
ends of tracks in Fig. [2] In some cases we either could not obtain 
initial convergence or could not properly maintain it until one of the 
criteria above was met. Results from such models were discarded. 

We performed control calculations on a single solar-mass 
star with different WIMP masses and annihilation cross-sections 
(Fig. |3). When capture and annihilation have equilibrated in a 
star, the WIMP luminosity effectively depends on the capture 
rate alone, which is independent of the annihilation cross-section 



(cf. Eqs. 2.14 2.15 and|2.19|, so we see that (a a v)o makes no dif- 



ference to the amount of energy generated. When equilibrium has 
not been achieved, this will not be the case. Such an effect can be 
seen in the slight upturn of the WIMP luminosity in the uppermost 
curve of the upper panel in Fig. [3] In this case the very high p x and 
very low (<T a i>}o significantly change the stellar structure before 
equilibrium has been reached, causing Lw to peak prior to equi- 
librium. The dependence of Lw on m x in Fig.|3]shows roughly an 
inverse square relationship, which is a result of using the full cap- 
ture expressions in the RSC. As can be seen from careful inspection 
of Eq. |2.19| for example, for small v and such as those used in 
the context of the early universe by e.g. |Iocco] (2008) and |Freese| 
|et al.| ( [2008) , the dependence disappears. 

In Figs. [4] and [5] we show evolutionary tracks in the HR and 
central equation-of-state diagrams of stars with different masses 
and WIMP luminosities. At low WIMP luminosities, the evolution 
is essentially normal. As WIMPs are allowed to provide more en- 
ergy, the negative heat capacity of a star causes it to expand and 
cool. The central temperature and density drop, nuclear burning re- 
duces and the star moves some distance back up the Hayashi track. 
The reduction in central temperatures and overall luminosities pro- 
vided by pp-chain and CNO-process hydrogen burning are illus- 
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Figure 4. Evolutionary tracks followed in the HR diagram by stars of various masses, when WIMPs provide different fractions of their total energy budgets. 
Filled, unlabelled circles indicate the starting points of tracks, whilst labelled ones give indicative ages during the evolution of 1.4Mq stars. Tracks have been 
halted when the star exhausts the supply of hydrogen in its core or reaches the current age of the universe. Stars with a greater luminosity contribution from 
WIMPs push further up the Hayashi track and spend longer there before returning to the main sequence. Stars which come to be entirely dominated by WIMP 
annihilation (bottom right) evolve quickly back up the Hayashi track and halt, holding their position in the HR diagram well beyond the age of the universe. 



trated in Fig. [6] These values are taken at the time £ a djust when a 
star has completed its initial reaction to the presence of WIMPs, 
which corresponds to the central temperature and density reaching 
their minima and the star arriving at the bottom-leftmost point of 
its travels in Fig. [5] At very high WIMP luminosities, the stellar 
core expands and cools drastically, moving stars a long way back 
along the pre-main sequence and effectively shutting down nuclear 
burning all together. Such an object becomes a fully-Hedged dark 
star, powered entirely and perpetually by WIMP annihilation. 

At intermediate WIMP luminosities, nuclear burning is sup- 
pressed rather than completely extinguished. Its continued contri- 
bution to nuclear processing slowly raises the core temperature and 
density once more, in turn increasing the rate of nuclear reactions 
and accelerating the process. The star burns hydrogen alongside 
WIMPs, and goes on to evolve through a hybrid WIMP-hydrogen 
main sequence. Such evolution can be best seen in the bottom-left 
panel of Fig. [4] Thanks to the energy input from WIMP annihi- 
lation, the time it takes such a star to consume its core hydrogen 
is lengthened, so its effective main-sequence lifetime is extended 
(Fig. 17). The increase in main-sequence lifetime is notable at all 
metallicities, but most prominent at low Z, essentially because nor- 
mal main-sequence lifetimes are shorter at lower metallicity. We 



did not see changes with metallicity in the central temperatures, 
pp-chain or CNO luminosities of the stars in our grid. 

We should point out here that in the extreme case of a very 
large WIMP luminosity, it is highly questionable whether a star 
would have ever reached the main sequence at all, or if it might 
have simply halted during its initial descent of the Hayashi track. 
The same might even be true of stars with intermediate WIMP lu- 
minosities, since it is not clear whether nuclear burning will win 
out over WIMP annihilation at exactly the same ages and capture 
rates when a star is evolved from the main sequence as when it is 
evolved from the pre-main sequence. This behaviour has been seen 
explicitly by |Iocco"eT al. (20081. We strongly suspect that the so- 
lutions we find for high WIMP luminosities are the same as those 
obtained with models begun from the pre-main sequence by |Iocco| 
|et al.|(2008) . We are currently investigating this question in detail. 

The greater influence of WIMP capture and annihilation upon 
low-mass stars is strongly apparent in Fig.|7J Even if WIMPs sup- 
ply only a tenth the energy of nuclear burning, the lifetime of a 
0.8 Mq star is increased by almost a billion years. With the same 
ratio of WIMP-to-nuclear burning, the lifetime of a 2.0 Mq star is 
unchanged. Considering that according to Fig. [2] roughly an order 
of magnitude more dark matter is required for a 2.0 Mq star to even 
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Figure 5. Evolutionary tracks followed in the central equation-of-state diagram by the stars of Fig. [4] Filled circles indicate the starting points of the tracks. 
Dashed lines show the location at which hydrogen burning becomes the dominant energy source, which defines the zero-age main sequence (ZAMS). To 
the bottom-left of this line, the core is too cool and diffuse to support the star by nuclear burning alone. Stars with a greater luminosity contribution from 
WIMPs push further into this region as they ascend the Hayashi track and their cores cool and expand, and thus take longer to recontract and return to the main 
sequence. The slight departure from smoothness apparent in some curves is simply due to finite temporal resolution of models. 



achieve the same WIMP-to-nuclear burning ratio as a 0.8 Mg star, 
lower stellar mass is clearly a highly favourable property in the ob- 
servational search for dark stars. 

In Fig. [8] we show the extent of convection at t = t a djust 
in stars of various masses, as the WIMP luminosity is increased. 
Because WIMP annihilation is far more concentrated at the cen- 
tre of a star than nuclear burning, stars with higher WIMP-to- 
nuclear burning ratios exhibit steeper radiative temperature gradi- 
ents in their cores. This produces convective cores of increasing 
size, as the height over which the temperature gradient is supera- 
diabatic increases. In parallel, the overall cooling and expansion 
of the star results in cooler surface layers, increasing the H con- 
centration and opacity and resulting in progressively deeper sur- 
face convection zones. At high enough WIMP luminosities the 
two zones meet and the star becomes fully convective. At lower 
metallicities, the promotion of convection is deferred until sig- 
nificantly higher WIMP luminosities, with the effect strongest in 
higher-mass stars. As an example, a 0.4 Mq star at Z — 0.02 be- 
comes fully convective at log 10 [Lw,max/£nuc(0)] ~ —0.1, and 
requires log 10 [L w ,max/£nuc(0)] w 0.3 at Z = 0.0003. A 1.4 M Q 
star on the other hand requires log, [L W .max /Lnuc(0)] « 0.8 at 
Z = 0.02, but log 10 [L w ,max/Lnuc(0)] > 2 at Z = 0.0003. 



We plot the dimensionless WIMP conductive effectiveness <£ 
(Eq. \2A6\ at t = tadjust f° r the full range of WIMP luminosities 
and stellar masses in Fig. [9] Due to the Knudsen-dependent and 
radial suppression factors (Eqs. |2.4l"| and |2.44] l, the contribution of 
WIMP conductive energy transport turns out to be small over most 
of the the parameter space. Eq. |2.46| is, however, a rather coarse 
measure of the significance of conductive transport by WIMPs; a 
more detailed comparison would investigate its role in the effec- 
tive opacity of nuclear matter in the star. At lower metallicities, the 
values of <£ become slightly larger, with the effect increasing with 
stellar mass. WIMP energy transport might therefore be worth in- 
cluding in future studies of PopIII dark star evolution. 

A number of numerical features complicate the interpretation 
of Figs.[6]j9] To be able to simulate a grid of ~3500 stars and pro- 
duce a manageable amount of output data, we chose only to save 
model data every tenth timestep. This meant the resolution avail- 
able for choosing the points at which to define Lw.max and t a djust 
was not as high as it could have been. The effect upon i a djuat was 
greater than on Lw.max, as it was compounded by the fact that 
tadjust is not as simply found as its definition would have one be- 
lieve. In cases where the equilibration time-scale is long, t a djust can 
be comparable to or longer than r cq . When this happens, £ a djust and 
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Figure 6. Central temperatures (top) and total luminosities provided by hy- 
drogen fusion via the pp-chain (middle) and CNO-process (bottom), as a 
function of the luminosity provided by WIMP burning. Central tempera- 
tures and hydrogen-burning luminosities are as measured at t a dj us t> the 
point at which a star has just completed its initial adjustment to the presence 
of WIMPs in its core. This corresponds to the time at which models' evo- 
lutionary paths have reached their bottom-leftmost points in Fig. [5] Rates 
of energy production from hydrogen burning are expressed as fractions of 
their initial values. 



r cq begin to lose meaning, as the adjustment alters the capture rate, 
which feeds back on the adjustment. This only occurs when C(t) 
and A(t) are of intermediate size, because even though r cq is at its 
longest when C(t) and A(t) are very small, very little adjustment 
is necessary in this case so tadjust is also very small. The upshot of 
all this is that a small amount of noise appears in Figs. [6] [8] and [5] 
For the sake of aesthetics, we recomputed individual evolutionary 
tracks in Figs.|4]and|5]with data saved every timestep. 

For 1.0 and 1.2 Mq stars, the reference CNO-process lumi- 
nosity extracted at t = in Fig. [6] was somewhat overestimated, 
due to the initial relaxation of the stellar models. The CNO process 
is only just present in 1 .0- 1 .2 Mq stars, and extremely temperature- 
sensitive, so is very likely to be significantly altered during numer- 
ical relaxation. The overestimation is the reason curves for 1.0 and 
1.2 Mq do not tend properly to zero at low WIMP luminosities in 
the lower panel of Fig. [6] 

Despite extensive prior testing, we found that our stopping cri- 
terion ii) was sometimes not quite stringent enough. Occasionally, 
models were halted which would have just managed to leave the 
main sequence in less than the age of the universe. We removed a 
small number of stars we suspected this of having influenced from 
the grid. As a result, the exact slopes of the steepest parts of the 
curves for higher masses in Fig.|7]are somewhat uncertain. 

Some noise also exists in the plots of Fig. [7] simply due to 
the finite temporal resolution of the models. Timesteps typically 
become longer once a star nears the end of the main sequence, so 
some temporal 'overshoot' can occur before criterion i) is triggered. 
As always, our choice of the internal timestep scaling was a com- 
promise between obtaining the smoothest results and being able to 
compute a reasonable number of models in a tractable timeframe 
(proper treatment of WIMP capture makes dark stellar evolution 
far more time-consuming than standard evolution). 

For a small window of WIMP densities, we also found that 
stars underwent seemingly random expansion and recontraction 
events during their evolution on the hybrid WIMP-hydrogen main 
sequence. Suggestively, these windows correspond approximately 
to the areas where we were unable to find solutions using our static 
code iFairb atrn et al.|20 08). The results of the static code suggest 
that this window can be thought of in the following way: for a par- 
ticular ambient density of WIMPs, if one tries to obtain a static 
solution, one may find that two solutions exist for a WIMP-burning 
star. The first solution corresponds to a star where the central tem- 
perature is rather close to that of a normal star of the same mass, 
such that the WIMPs are spread over a larger volume according to 
equation |2.25| Because the WIMPs are spread over a larger vol- 
ume, their annihilation is less centralised and the spatial distribu- 
tion of their energy release into the star is closer to that of normal 
nuclear burning. The core temperature of the star is therefore not 
reduced very much and the solution is self consistent. The second 
solution occurs when the WIMPs are localised in a small region in 
the centre of the star because the central temperature is low. The 
low central temperature is what would be expected if energy were 
input into a small region at the centre of the star, so this solution is 
also self-consistent. This situation arises primarily for larger mass 
stars, partially due to the extreme temperature-dependence of the 
CNO process which dominates nuclear burning in such stars. 

The existence of two solutions with the static code suggests 
that we might expect evolution in this region of parameter space to 
be unpredictable when looked at with a time-dependent code. For 
example, stars might exhibit genuine periodic or chaotic variabil- 
ity, or could become numerically unstable. By reducing the internal 
timestep scaling and employing the reconvergence mode, we es- 
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Figure 7. Main-sequence lifetimes of stars with differing masses, metallicities and WIMP luminosities. Dashed lines indicate the present age of the universe. 
Stars with a greater fraction of their energy provided by WIMPs consume hydrogen more slowly, so spend more time on the main sequence. For a given 
WIMP-to-nuclear burning ratio, the lifetimes of stars with lower masses and metallicities are more affected than their more massive, metal-rich counterparts. 



tablished that the apparently random expansions and contractions 
are numerical artefacts caused by insufficient temporal resolution; 
given sufficiently small timesteps, the time-dependent code seems 
to follow a path intermediate to the two solutions appearing in the 
basic static code. The borderline stability of the stellar structure in 
this region when treated with a code which assumes hydrostatic 
equilibrium (as DarkStars does) suggests that such stars might ex- 
hibit some true physical variability after all, but due to dynamical 
effects only captureable with a full hydrodynamic code. If there is 
interesting variable behaviour in this region then, it probably occurs 
on a timescale smaller than can be resolved by DarkStars. 

The evolution of such stars over an entire lifetime seems 
largely unaffected by the excursions, so we are confident that the 
overall results of the grid still accurately reflect the general proper- 
ties of main sequence dark stars. However, the excursions do fur- 
ther complicate the task of automatically choosing t a dj us t. adding 
further noise to Figs.|6][8]and|9] 



4 THE GALACTIC DARK MATTER HALO 

In order to simulate the accretion of dark matter by stars at the 
Galactic centre, one should understand its distribution in the Milky 



Way. In particular, the density and velocity distribution of dark mat- 
ter both play a role in determining the capture rate. 



4.1 Density 

Density profiles of dark matter halos have been a topic of compu- 
tational study for over a decade l Navarro et al.| 1996} |Moore et al.| 



[T999"l|Navarro et al.|2004"l|Diemand et al.|2007) . As TV-body sim- 
ulations have been run on computers of ever-increasing speed, a 
standard lore for the expected distribution of dark matter in halos 
of all sizes has developed. Two conclusions that seem to be uni- 
versally accepted are that dark matter is denser in the centre of 
halo simulations, and that the logarithmic gradient of the density 
7 = dlnp/dlnr is more negative in the outer parts of simulated 
halos than the inner parts. A popular parametrisation is the 'NFW 
profile,' which interpolates smoothly between asymptotic values of 
7 = — 1 in the inner regions of the halo and 7 = —3 in the outer 
regions ( N avarro et al.|1996} . |Navarro et al.| l |2004) have since sug- 
gested that a better profile would in fact be one where 7 varies 
smoothly with radius, and does not asymptote to any particular 
value at large nor small radii; such a profile has been advocated 
by various authors since the 1960s (see M erritt et al.|2006] for an 
historical account). 

Most simulations only include dark matter, since considering 
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Figure 8. The extent of convection at t = t a dj us t m stars of different masses as the energy from WIMP annihilation in their cores is increased. Shaded 
areas indicate regions in which the stellar energy transport is convective; elsewhere, transport is radiative. Stars develop and extend their convective cores 
and envelopes as the WIMP luminosity is increased, eventually becoming fully convective at high values of max - Less massive stars require a smaller 
ratio of annihilation to fusion energy to alter their convective properties than heavier stars. Plots in different panels extend over slightly different ranges of the 
WIMP-to-nuclear burning ratio, according to which models we were able to converge. Not shown is the dependence upon metallicity: at lower metallicities, 
the onset of convection is deferred to higher WIMP luminosities, with effects greatest in higher-mass stars. 



collisionless particles does not require the complicated hydrody- 
namics necessary to model baryons. The presence of baryons is ex- 
pected to change the distribution of dark matter: since the baryons 
are able to lose energy end sink into the middle of a galaxy, they 
create a potential well which subsequently pulls the dark matter 
into the central region. This phenomenon of adiabatic contraction 
was first predicted by Zel'dovich et al.| < j 1980| >, and formalised by 
|Blume nthal et al. ( 1986 1. More recently, it has been realised that the 
non-circular nature of typical dark matter orbits reduces this effect, 
but does not remove it ( |Gnedin et al. 2004). In order to calculate 
the expected density profile of dark matter in a given galaxy, it is 
therefore necessary to take a dark matter halo from an TV-body sim- 
ulation of an appropriately-sized galaxy, then adiabatically contract 
it using the galaxy's observed baryonic profile. This should be done 
taking into account the non-circular nature of the dark matter orbits 
( |Gnedin et al.|2004[|Gustafsson et al.|2006) . 

This prescription gives a more realistic estimate of the ex- 
pected dark matter density in the central regions of a real galaxy 
than the results of a collisionless TV-body simulation (for a compar- 
ison of results obtained using this procedure with pure dark matter 
halos and those from simulations which also include baryons, see 
|Gustafsson et al.||2006[ >. The typical effect of such a contraction 
is to draw dark matter deeper into the central part of the galaxy, 
changing the inner slope of the density profile. 

As one approaches the centre of a galaxy like the Milky Way 



from a large distance, the gravitational potential is first dominated 
by the diffuse dark matter halo. Approaching the central bulge, the 
gravitational potential of the concentrated baryonic mass becomes 
more important. If the current understanding of the dark matter dis- 
tribution in the Milky Way is correct, the changeover occurs at 
a radius of the same order of magnitude as the solar position. In 
the centremost regions, the supermassive black hole determines the 
gravitational dynamics. The density of dark matter is thought to rise 
continuously towards the centre of the galaxy, but at radii much less 
than the solar position, its gravitational influence is always dwarfed 
by that of baryons or the central black hole. 

In the central parsec, where the black hole starts to dominate 
the gravitational potential, the dark matter profile depends upon a 
number of factors. One is simply the density of dark matter at larger 
radii, which forms an initial condition for the central density pro- 
file. If the black hole forms in situ, it may create a miniature adia- 
batic contraction of the local dark matter profile, leading to a cen- 
tral spike where the density gradient is steeper than in the rest of the 
galaxy. Immediately after its formation, the density of dark matter 
in such a spike can be extremely high (Gondolo & Silk 1999). The 
spike is expected to diffuse away over time due to dark matter self- 
annihilation, loss of dark matter as it falls into the black hole and 
heating of the dark matter by gravitational interactions with stars. 

These different effects can be incorporated into a diffusion 
equation that gives rise to a final prediction for the density of 
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Figure 9. The significance of conductive energy transport by WIMPs, as 
measured by the dimensionless WIMP conductive effectiveness (Eq. |2.46) . 
A value of on the y-axis roughly indicates that conductive energy transport 
by WIMPs is as important as all other energy sources in the star combined. 
The importance of conduction by WIMPs is significantly less than that of 
actual energy sources. Not shown is the fact that £ becomes slightly larger 
at lower metallicity, with the effect increasing with stellar mass. This sug- 
gests that conductive energy transport should probably be taken into account 
when simulating massive metal-free stars. 



Table 1. Parameters for the density profiles defined by Eqs. |4.1| which are 
approximations to the profiles presented by Bertone & Merritt[j2005). 



Profile 


p x (100pc) 


71 


72 


rout 


Tin 


NFW+spike 
AC+spike 


25GeV cm" 3 
360 GeVcm" 3 


1 

1.5 


1.85 
1.82 


7 ■ 10 4 R BIl 

7 ■ 10 4 R BH 


10 R BH 

IORbh 



dark matter near the central black hole (Bertone & Merritt 20051. 
In this work we will consider two dark matter density profiles, 
both approximations to the profiles presented by Berto ne & Merritt| 
( |2005|>. These appro ximations correspond to profiles B and C used 
by |Ber gstrom et ai~| ( |2006fr . The 'NFW+spike' profile is a standard 
NFW o = l profile with a central spike which has diffused away 
over time, considerably reducing its density. The 'AC+spike' pro- 
file has undergone adiabatic contraction on galactic scales due to 
the presence of baryons, and also has a central spike which was al- 
lowed to diffuse away over time. Both profiles can be parametrised 
by the expressions 

/ a x (r) = Px (100pc)(i^) 71 r>r out 

Px(r) = Px(r a ut) ( I v i ) 72 r out >r>r ln (4.1) 

Px( r ) = Px( r in) r ln >r, 

where parameters are listed in Table [T] Note that following adi- 
abatic contraction, the smoothly varying profile advocated by 
|Navarro et al.| ( |2004} can become almost as steep in the central 
region as an equivalently-contracted NFW profile, depending upon 
the angular momentum of the dark matter particles. 



4.2 Velocities 

Having obtained some estimates of the possible densities of dark 
matter at the Galactic centre, we must also think about its velocity 



distribution, which has a strong bearing upon the number of parti- 
cles captured by stars. 

Various estimates of the velocity distribution exist in the lit- 
erature. For direct detection experiments such as CDMS, Xenon 
and COUPP to be able to easily compare results, the dark mat- 
ter halo is typically assumed to be the isothermal RSC, with a 
radius-independent Keplerian velocity. In this case, the velocity 
dispersion is set by the Keplerian velocity at the solar position 
(v — ^/3/2dq = 270kms _1 ), and acts as the width for a Gaus- 
sian distribution of velocities which is identical at every position. 

As already mentioned, the highest resolution V-body simula- 
tions do not predict an isothermal halo, but rather one where the 
logarithmic density gradient close to the centre of the galaxy is 
much less than —2. The Keplerian velocity in such a halo would 
therefore decrease to zero at the core, which would increase the 
rate at which dark matter would be accreted by stars. In a real 
galaxy however, the dark matter is a subdominant component at 
these small Galactic radii, and the presence of stars and the central 
black hole increases the gravitational potential there. 

Our default assumption is that the velocity distribution of dark 
matter is isotropic, spherically symmetric, Gaussian and has a dis- 
persion set by the Keplerian velocity in the solar vicinity. None of 
these assumptions are strictly correct, as we shall discuss shortly. 
Our simplest attempt to improve the realism of the velocity dis- 
tribution is to exclude velocities above the local Galactic escape 
velocity, as WIMPs with such velocities would presumably already 
have left the Galaxy some time earlier. We therefore truncate the 
velocity distribution at the local escape velocity (in the rest frame 
of the Galaxy), which terminates the Maxwell tail of the distribu- 
tion. Given a generic WIMP velocity distribution go [u], seen in the 
rest frame of the galaxy, the equivalent truncated distribution will 
be 



ff ga i,o(M) 



Px go{u)S{u sal - u) 
m x / °° go{u')0(u gal ~ u') du' 



(4.2) 



where w ga i is the local escape velocity, and the integral ensures 
the new distribution remains correctly normalised. For a star at rest 
with respect to the Galactic halo, this then brings the capture rate 
(Eq.|2.2} to the form 
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The normalisation factor D is 
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Working from Eq. |2.12| in the case of an isothermal velocity distri- 
bution this becomes 
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(4.5) 



To find the capture rate in the frame of a star moving relative 
the the Galactic rest frame, we need to transform g ga i,o (if) to some 
equivalent g ga i,*(w) via an appropriate Galilean transform, in anal- 
ogy with the step from Eq. |2.12| to Eq. |2.13| then consider what 
iigai becomes in the frame of the star. The maximum velocity any 
WIMP from a distribution cut off at u ga i can have in the galactic 
frame is obviously u ga i. In the frame of the star though, WIMPs 
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coming from e.g. the direction in which the star is moving through 
the halo can appear with much greater speed than those coming 
from 'behind' the star. If an incoming WIMP has speed uo and ve- 
locity in the direction of the unit vector ew in the galactic frame, 
we see that its velocity in the star's frame is 

= Uo€w — V* 

=> u = = v (tmew - v*) 2 

U < \J (ll ga iew — V*) 2 = (llg al + v'l + 2u gal V ic COS If) 3 . 

So in the star's frame, we have 



Mgal 



Wgal,*(<y3) = (^gal + V* + 2u ga .\V ir COS if) ■ 



(4.6) 



where if is the angle in the galactic frame between the motions 
of the WIMP and star, such that f — corresponds to a head- 
on collision. This then poses something of a problem, as in the 
analogous case of the isothermal distribution, if had to be implicitly 
integrated over in the first place to obtain the integrand (Eq. |2.13} 
for which this new cut-off velocity is the limit. The solution to this 
is either to develop some sort of approximate averaging scheme, 
or to just do the integral over if explicitly, after the integral over 
u. In this case, the expression for capture of a truncated isothermal 
distribution of WIMPs by a moving star is 



C(t) = 



P^D- 1 I r 2 
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x wQ v i (w) exp 
+ 2uVi, cos tp) J du d(cos ip) dr, 

where 

u cu t,i(f,r,t) = min[M ga i,*(y)) I w Ina x,i(r, t)]. 



(4.7) 



(4.8) 



As in the non-truncated case, the velocity integral in Eq. |4.7| 
can be performed analytically, yielding the truncated analogues of 
Eqs. |2~T4}(2~T9l 

C(t) = 47TD- 1 / r 2 Y] [Mi(« eu t, i ((p,r,t)) 
Jo J-l i 

- Mi (0)] d(cos^)dr, (4.9) 
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for heavier elements, and 
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for hydrogen. This gives the machinery necessary to calculate the 
capture rate for a truncated isothermal (Gaussian) velocity distribu- 
tion. 

The fact that we do not expect the velocity distribution to be 
Gaussian in reality is related to the fact that one only expects a 
Gaussian distribution in the limit of an extensive distribution of 
particles, such as an ideal gas. For particles coupled to a long-range 
potential such as gravity, this is not the case. It has been suggested 
in the literature that the Tsallis distribution, one designed specifi- 
cally to model the departure from extensivity, is a better fit to the 
data than a Gaussian (Hanse net al.|2006) . 

Furthermore, there are good reasons to believe that the radial 
distribution of dark matter will have a different width to the tan- 
gential distribution, since the orbits of dark matter particles are far 
from circular. For a star on a non-circular orbit around the centre of 
the Galaxy, this could have important consequences. 

Finally, it is interesting to test the truth of the assumption that 
the velocity dispersion is fully determined by the Keplerian veloc- 
ity at the solar position. The relationship between the Keplerian ve- 
locity, the velocity distribution and its anisotropy depends upon the 
shape of the potential well of the galaxy; the anisotropy in the ve- 
locity distribution can be obtained from the Jeans equations ( |Fair-| 
|bairn & Schwetz|2009) . 

In order to quantify the departures from the isothermal halo 
model, we have examined data from the Via Lactea simulation 
( |Diemand et al.|20 07 1. This TV-body simulation contains more than 
2 x 10 8 dark matter particles, and is one of the largest simulations 
of a Milky Way-size dark matter halo to date. As suspected, the re- 
sults do indeed show that all four of the simplifying assumptions 
involved in the isothermal halo model (isotropy, spherical symme- 
try, Gaussianity and a dispersion proportional to the Keplerian ve- 
locity at the solar position) are essentially incorrect. To obtain a 
new velocity distribution from the data, we looked at the velocities 
of particles at different radii. We attempted to fit the distributions 
with the Tsallis profile of Hansen et al. (2006 ), but found that al- 
though this does provide a better fit than a Gaussian, the following 
one-dimensional distribution is an even better fit: 



hm(iii) = exp < - 



1 fu^ 2 
2\<Ji 



(4.13) 



where i G {r, 0, <f>} and Ui is the velocity in the zth direction in 
Galactic coordinates. No normalisation pref actor has been included 
here. The parameter a measures the departure from Gaussianity 
and is distinctly different for the one-dimensional velocity distribu- 
tions in the radial and angular directions. The values of a for the 
radial distribution (a r ) and a composite tangential distribution (ot, 



where u T 



: Ug + w|) can be found as a function of radius in 



Fair- 



|bairn & Schwetz ( 2009 ). The ratio of the velocity dispersion to the 
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local Keplerian velocity (uk cp , or more precisely, the square root 
of the potential since the halo is not completely spherically sym- 
metric) is also a function of radius and can be found in Fairbairn & 
|Schwetz|j2009) . 

These distributions become less Gaussian as one approaches 
the Galactic centre, with a <C 1 in the region we are most inter- 
ested in. One should be careful not to take this result overly seri- 
ously though; not because the physics of the simulation is in doubt 
(above the ~0.1 kpc resolution scale), but rather because the sim- 
ulation only includes dark matter. The velocity dispersion at the 
centre of an NFW profile goes to zero, whereas in a real galaxy, the 
presence of baryons and the central black hole would be expected 
to change the velocity distribution quite significantly. We do not 
make any strong claims as to the fitness of this non-Gaussian ve- 
locity distribution for modelling the very centre of the Galaxy; we 
employ it more with the goal of determining what degree of uncer- 
tainty exists in our capture results due to the velocity distribution. 

In order to reasonably calculate the effect of the non-Gaussian 
distributions upon capture rates, we require a composite distribu- 
tion of total velocity magnitudes in three dimensions. Since the fit- 
ted values of a and ct/vkcp become almost isotropic at low Galac- 
tic radii, for this purpose we can set them to the same values in 
the radial and angular directions. We choose these values as the 
means of the fitted values in each spatial direction, for the small- 
est Galactic radius at which we fit the velocity distribution (1 kpc). 
This gives a — 0.35 and <t/«k cp = 0.05. At the Galactic centre, 
VKep is dominated by the black hole, so 

a = 56.8kms- 1 x j^ ^) 172 . (4.14) 

To obtain a three-dimensional distribution purely as a function of 
the velocity magnitude 



y^ 2 + u 2 + 



(4.15) 



we take the product of the 3 one-dimensional distributions, convert 
to local spherical polar coordinates and integrate over the angular 
directions. This gives 



hsD(u r , ue, u^) = iiru exp 
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When a — 1, this clearly produces a Gaussian distribution of ve- 
locity magnitudes. For a / 1, we can express u 2 ,, Ug and w| in 
terms of one another and it 2 , then expand to second-order in a bi- 
nomial series to produce 
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(u 



2\a 2a 
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(«; +Uj), 



(4.17) 



where i and j are different members of the set {r, 9, </>}. This gives 



We now have a three-dimensional velocity distribution as a function 
of the local Keplerian velocity, which can be inserted into Eq. |2.2| 
to obtain the capture rate. Putting Eq. |4.19| into Eq. |4.3| instead, one 
obtains the capture rate from an equivalent truncated distribution, 
which is also a function of the local escape velocity. 



5 THE GALACTIC POTENTIAL 

To calculate dark matter capture rates, it is extremely important to 
know the stellar velocity v* through the dark matter halo. In order 
to correctly truncate the velocity distribution of the dark matter, one 
also needs to know the local Galactic escape velocity u max . 

The orbital velocities are rather simply obtained. Within a few 
tenths of a parsec, the Galactic potential is dominated by the central 
black hole. All the elliptical orbits we will consider lie within one 
fiftieth of a parsec, so they can be treated as exactly Keplerian about 
a point-mass black hole. 

Calculating the escape velocities is more arduous, as we need 
to integrate the potential experienced by a test particle exiting the 
Galaxy. It is important to consider not only dark matter but also 
the presence of baryons, which dominate the potential from around 
0.5 pc to several kpc. To model the baryon density of the Milky Way 
we use the same prescription as Gust afsson et al.| ( [2006| l, assuming 
a central bulge of stars with density p oc r 1 e~ r/A . We assume a 
thin disc of matter with surface density 
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cM dh 



2tt (r 2 + c 2 )3 



(5.1) 



We choose the free parameters to match observations of the Milky 
Way: 7 = 1.85, A = lkpc, c = 5 kpc and the tot al disc and 
bulge masses M d i SC oo = 5M bu i ge = 6.5 x 1O 1O M jKent et al. 
p9T] [aiaol[T996l |Detmen & Bumeyl[T998l |Klypin et al.||2002| >. 
We assume that the extent of the disc is 15 kpc. We use an NFW 
profile with a scale radius of 20 kpc for the dark matter, normalised 
to 0.3GeVcm -3 at the location of the solar system, and find the 
local escape velocity by integrating the energy loss along a radial 
path exiting the Galaxy. 



6 IMPACTS ON STARS AT THE GALACTIC CENTRE 

Armed with detailed estimates of the stellar orbits, local escape 
velocity, density of dark matter and its velocity distribution at the 
Galactic centre, we can now realistically evaluate the potential 
impact of WIMPs upon stellar evolution there. We ran two fur- 
ther grids of evolutionary models, both at Z — 0.02 and over 
0.3 M Q ^ ^ 1.5 M . We computed models in both grids 
using the NFW+spike and AC+spike profiles from Sect. [4] 
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Normalising over u £ [0, oo) gives the final three-dimensional 'TV- 
body' velocity distribution 
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Whilst this expression breaks down at third order in the binomial 
expansion, we expect it to be a reasonable approximation given the 
level of uncertainty in choosing realistic values of a and a/vKe P - 



6.1 Circular orbits 

The first grid covered single stars on circular orbits, with orbital 
radii extending from 10 pc to 10 _6 pc. For this grid we used the 
standard, non-truncated version of the isothermal velocity distribu- 
tion (Eq. |2.13[ >. Results for models computed with the AC+spike 
density profile are shown in Fig. [TU] As orbits are made smaller, 
capture rates rise because stars encounter higher densities of dark 
matter. This effect is balanced by a reduction in capture caused by 
stars' increasing circular velocities as they orbit closer to the black 
hole. The capture rate and resultant WIMP luminosity peaks at an 
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Table 2. Orbits considered in Sect. |6.2| along which we evolved stars in 
Figs.[Tl][l3] 



-4 -3 -2 
Orbital radius, log 10 { 



Figure 10. WIMP luminosities achieved by stars orbiting circularly about 
the central black hole. The dark matter velocity distribution is isothermal 
with dispersion v = 270 km s - , and the density profile follows the adia- 
batically contracted profile (AC+spike). The impact of instead using a ve- 
locity distribution truncated at the local Galactic escape velocity is shown 
with a dashed curve. Capture is maximised at a radius of approximately 
0.3 pc, but no circular orbit produces capture rates high enough to trans- 
late into a WIMP luminosity which can produce any significant changes in 
stellar evolution. WIMP luminosities produced with the alternative density 
profile (NFW+spike) are even lower, so are not shown. 



orbital radius of ~0.3pc. Inwards of this the velocity effect domi- 
nates and capture is highly suppressed. The constant WIMP lumi- 
nosities at very small radii in Fig. [10] are entirely due to the initial 
populations of WIMPs that the models were started with. 

We evolved some supplementary models (dashed curve in 
Fig.| 1 0j> wit h the truncated isothermal distribution (Eq. |4.2| applied 
to Eq. 2. 12) , to see if capture might be boosted to interesting lev- 
els by removing unphysical WIMP velocities. On the contrary, the 
truncation of the isothermal WIMP velocity distribution caused a 
strong reduction in capture rates. Stars moving as quickly as those 
on circular orbits near a black hole must capture predominantly 
from the Maxwell tail of the isothermal distribution, so truncation 
denies them many of their best capture candidates. The opposite is 
true of a star in the RSC, which captures from the centre of the dis- 
tribution and benefits (slightly) in capture rate if the distribution is 
truncated and renormalised. 

If it follows a circular orbit in an isothermal WIMP halo, even 
the lowest-mass single star, placed at the optimal radius in the most 
optimistic density profile, cannot accrete enough WIMPs to bring 
annihilation luminosity to within two orders of magnitude of its 
nuclear luminosity. We do not show WIMP luminosities resulting 
from the NFW+spike density profile, as they are even less interest- 
ing due to the profile's lower central density. 



6.2 Elliptical orbits 

The primary assumptions of the previous grid were that stars al- 
ways follow circular orbits, and that the halo is isothermal. We 
know these to be untrue in reality, so in the second grid we con- 
sidered the effects of elliptical orbits and a non-Gaussian velocity 
distribution. 

This grid consisted of single stars on orbits with various ellip- 
ticities, within three classes: orbits with periods P — 10 yr, orbits 
with P — 50 yr and orbits where the maximum star-black hole 
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separation was 0.01 pc. The galactocentric distances at periapsis 
and apoapsis (r m in and r max respectively) of each of the orbits we 
considered for this grid are given in Table |6. 2] The maximum ec- 
centricity in each class was chosen so as to ensure that stars did not 
come within five Scwharzschild radii of the centre of the black hole, 
ensuring that relativistic corrections to the orbit were not critical. 

The early evolution of one of the stars from the grid is shown 
in the left panel of Fig. [TT] The initial flat sections of curves are 
where the star was held at a constant galactocentric radius for the 
first 3 timesteps to allow the model to properly relax. Capture and 
annihilation occur in punctuated stages, clearly correlated with the 
orbital period. Strikingly, the times of greatest capture are in fact 
when the star is farthest from the centre of the Galaxy, at apoapsis. 
This is because it has had a chance to slow down relative to the DM 
halo, and achieve a significant capture rate for a time before plum- 
meting back down towards the black hole. By the time it reaches 
periapsis, the star is moving so quickly that capture is essentially 
zero, regardless of how high the dark matter density is. 

Because following the dark evolution of a star on such an or- 
bit is extremely time-consuming, we evolved the models in this grid 
for just 5 full orbits each, then calculated the average capture rates 
achieved over this time. We explicitly assume that given identical 
initial mean capture rates, the long-term evolution of a star on a 
short-period elliptical orbit would be the same as the evolution of 
one which evolves on a circular orbit. That is, we assert that be- 
cause of the one-to-one mapping between capture rate, WIMP lu- 
minosity and evolutionary effects discussed in Sect. [3] the evolution 
of stars on elliptical orbits in the Galactic centre can be predicted 
by assigning them an equivalent RSC star from the grid of models 
we presented earlier. Recall now our assumption in Sect. |2.2| that 
WIMPs instantaneously thermalise with the stellar material. This is 
almost certainly not a good approximation on the timescale of just 
5 orbits. We therefore probably overestimate annihilation rates and 
how closely they track capture during such stars' early years, since 
annihilation takes longer to catch-up with capture when instanta- 
neous thermalisation is not assumed. As our primary goal is to pre- 
dict the long-term evolution on the basis of the initial capture rate, 
we don't expect this approximation to have a large impact here. We 
plan to directly test all these assumptions in later papers, using sin- 
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Figure 11. Evolution of a 1 Mq, Z = 0.02 star on a highly elliptical orbit close to the central back hole, followed for 5 orbits in the beginning of its lifetime 
(left) and at an age of half a billion years (right). The orbit is that listed in the final line of Table |6.2] Dark matter velocities follow an isothermal distribution 
with v = 270 kms -1 , truncated at the local escape velocity, and densities follow the AC+spike profile. WIMP capture occurs exclusively around apoapsis 
despite this being the point at which the ambient WIMP density is lowest, thanks to the stars very low orbital velocity. Early in the evolution, before capture 
and annihilation have equilibrated and the WIMP population has stabilised, the total population, annihilation rate and resultant WIMP luminosity undergo 
punctuated increases each time the star goes through a period of capture. In the evolved star on the right, the equilibrium population of WIMPs provides a 
buffer against the transient nature of capture, and the evolution is essentially smooth. Because explicitly following the evolution on an elliptical orbit is highly 
time-consuming, between these two plots the star was allowed to evolve on an artificial circular orbit with the same initial mean capture rate as exhibited in 
the left panel. Because of the finite temporal resolution, some of the peaks at periapsis are not properly resolved; because capture here is effectively zero, this 
has no effect upon the evolution. 



gle simulations with very long runtimes and explicitly simulating 
the thermalisation process. 

Under these assumptions, we evolved the star of Fig. 1 1 1 1 for 
a further half a billion years with a constant equivalent RSC den- 
sity of 3 x 10 9 GeV cm" 3 , to let it fully adjust to the effects of its 
captured WIMPs. We then put it back on the same elliptical orbit, 
where it was allowed to evolve for a further 5 orbits (right panel 
of Fig.fTT). As an evolved dark star, it continues to capture in the 
same punctuated manner as during its early years. Now that the star 
has built up an equilibrium population of WIMPs though, its total 
population and annihilation rate are essentially immune to the tran- 
sient nature of capture. The small decrease in the total population 
and annihilation rate over the lOOyr shown here is just due to a 
slight mismatch between the chosen equivalent RSC density and 
the actual mean capture rate on the elliptical orbit. 

In an attempt to prevent very rapid changes in the capture rate, 
we limited timesteps to those which would prevent the ratio ^ 
shifting by more than 30 per cent in a single step. Because the 
evolution code has trouble with convergence when timesteps are 
reduced too far, we also had to impose a lower limit to timesteps 
of 0.5-1.7 yr in order to prevent the above criterion breaking con- 
vergence as stars passed through periapsis. Luckily, periapsis turns 



out to be the least important part of the orbit for the calculation of 
the mean capture rate and the actual evolution, so the only conse- 
quence of this is aesthetic: some of the peaks and troughs in Fig.|l 1| 
are not fully resolved. 

Fig. [12] shows the mean capture rates achieved by all stars in 
this grid of models, assuming a truncated isothermal velocity distri- 
bution. Thanks to the additional capture window opened at apoap- 
sis, stars following elliptical orbits have their capture rates boosted 
by up to 20 orders of magnitude beyond what they would have 
achieved on the equivalent circular orbit. The more elliptical the or- 
bit, the slower the star is moving at apoapsis, so the more dark mat- 
ter it is able to capture. Whilst stars in the two orbital classes with 
constant periods reach apoapsis further from the black hole with 
increasing orbital ellipticity, the reduction in capture rate caused by 
the resulting decrease in dark matter density at apoapsis is com- 
pletely outweighed by the increase in capture brought on by the 
reduced star- WIMP relative velocities. The ellipticity boost is most 
marked for shorter orbital periods, as stars on shorter-period orbits 
have the most to lose by following circular orbits (due to their very 
high circular velocities), yet the most to gain by following elliptical 
orbits (because they sample regions of higher dark matter density). 

Referring to Fig. [2] the stars of Fig. [12] evolved in the 
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Figure 12. Mean capture rates achieved by stars on elliptical orbits with 
periods of 50 years (top) and 10 years (middle), as well as orbits where the 
star-black hole separation is 0.01 pc at apoapsis (bottom). Dark matter ve- 
locities follow an isothermal distribution with dispersion v = 270 km s , 
truncated at the local Galactic escape velocity. Capture rates are boosted by 
up to 20 orders of magnitude when stars follow elliptical rather than circular 
orbits. 



Figure 13. Mean capture rates achieved by a 1 Mq star on elliptical orbits 
through dark matter halos with different velocity distributions. All halos fol- 
low the AC+spike density profile. The TV-body velocity distribution results 
in globally higher capture rates than the standard isothermal distribution. 
Truncating the velocity distribution at the local Galactic escape velocity de- 
creases capture rates from the isothermal distribution, but has no impact 
upon capture from the TV-body distribution. 
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Figure 14. WIMP luminosities achieved by stars on orbits with 10-year pe- 
riods around the Galactic centre. Annihilation can provide up to 100 times 
the power of nuclear fusion in stars on realistic orbits. If the Galactic halo 
has undergone adiabatic contraction (AC+spike), annihilation rivals nuclear 
fusion in stars on any orbit with an eccentricity greater than about e = 0.9, 
for all masses less than or equal to 1 .5 Mq . If not (NFW+spike), stars of a 
solar mass or less approach break-even between fusion and annihilation en- 
ergy on orbits with e > 0.99. These curves have been obtained by applying 
the boosts seen in Fig. |13| to the capture rates of Fig. [12] and then interpo- 
lating within the results shown in Fig. [2] to obtain the resulting WIMP-to- 
nuclear burning ratios. The arrow indicates that the 1 Mq, AC+spike curve 
is expected to continue in this direction, but there is no reliable way to con- 
vert capture rates to WIMP luminosities in this region because the capture 
rates and implied WIMP luminosities are beyond the range of convergence 
of the benchmark models in Sect. [3] In the case of the least interesting orbits 
(circular orbits in an NFW+spike density profile), capture rates and WIMP 
luminosities are in fact below the limits of the grid in Sect. [3] so WIMP lu- 
minosities are extrapolations based on the assumption that the scaling with 
ambient WIMP density is the same at e = as at e = 0.5. 
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AC+spike density profile can all achieve break-even between 
WIMP and nuclear luminosity if their orbits have a period of 10 yr 
and e > 0.99. A 0.6 Mq star can achieve this goal even on a 50- 
year orbit, with ellipticity as low as e = 0.9, whilst a 1 Mq star 
can do the same for e > 0.99. With the NFW+spike density profile 
and a truncated isothermal halo, stars will generally never achieve 
the same level of energy output from WIMP annihilation as nuclear 
burning (at least not if m x ^ 100 GeV; cf. Fig. [3}. 

In Fig. [T3] we show the impacts upon capture of going be- 
yond the isothermal halo approximation. Here we illustrate the cap- 
ture rates achieved by a 1 Mq star on 10-year orbits through both 
truncated and non-truncated versions of the A-body velocity dis- 
tribution, and compare with the corresponding capture rates from 
isothermal distributions. Regardless of the orbital ellipticity, cap- 
ture is at least a factor of 3 (i.e. half an order of magnitude) higher 
from the standard version of the A-body distribution than from 
the standard isothermal distribution. This difference blows signifi- 
cantly when the distributions are truncated at the local escape ve- 
locities and renormalised; whereas the truncation reduces capture 
from the isothermal distribution (particularly at low ellipticities), it 
leaves capture from the A-body halo entirely unaffected. 

The truncated A-body distribution boosts capture rates at high 
ellipticities by a factor of 3-5 over the truncated isothermal distri- 
bution, significantly increasing the range of orbits over which we 
might expect to see dark stars. In Fig.[l4]we estimate the long-term 
behaviour of stars on the 10-year orbits by combining the capture 
rates seen in Fig.[l2] the boosts seen in Fig.[T3]and the results from 
our benchmark simulations in Sec. [3] Here we have made the ap- 
proximation that boosts do not depend upon the stellar mass or halo 
density profile; that is, we adjusted the capture rates of all stars by 
the difference between the truncated A-body and truncated isother- 
mal curves in Fig. [13] depending only upon the orbits which stars 
followed. We then interpolated within the data of Fig. [2] to obtain 
WIMP luminosities from the adjusted capture rates, which included 
interpolating further amongst the curves of Fig.|2]to obtain data for 
1.5 M Q stars. 

Fig. [14] shows that the effects of WIMPs can be drastic, with 
the energy produced by WIMP annihilation outstripping that of nu- 
clear fusion by up to a factor of 100. Indeed, some of the adjusted 
capture rates imply a WIMP luminosity even greater than that of 
any star we were able to reliably evolve in the grid of benchmark 
models (indicated by an arrow pointing in the direction one would 
expect the curve to continue in). We now see from Fig.[l4]that even 
the NFW+spike density profile can actually produce stars where 
annihilation comes close to breaking even with nuclear burning 
(P w 10 yr, e > 0.99 and A/* < 1 M ). With the AC+spike pro- 
file, the same is true of stars with masses of up to 1.5 Mq, following 
orbits with eccentricities e > 0.9. If one were to perform the same 
conversion on capture rates achieved by stars on other orbits, pos- 
tulating similar boosts as seen on 10-year orbits and assuming an 
AC+spike density profile, stars of M* < 1 Mq would also achieve 
break-even for e > 0.9 with orbital periods of up to 50 years. 

Finally, we point out that the magnitude of these results de- 
pends upon the chosen WIMP mass (Fig. [5}. As such, one expects 
the final WIMP luminosities shown in Fig. [14] to be even further 
boosted for WIMPs lighter than 100 GeV, but suppressed for higher 
masses. The dependence of the capture rate and final WIMP-to- 
nuclear burning ratio upon the WIMP mass becomes weaker for 
smaller velocity dispersions, at least in an isothermal halo. Whilst 
we have not explicitly investigated how this picture changes in a 
halo with a non-Gaussian velocity distribution, we would expect 
similar behaviour. With velocity dispersions given by Eq. |4.14| we 



therefore expect WIMP luminosities in the A-body distribution to 
have some dependence upon the WIMP mass, though not so pro- 
nounced as seen in Fig. [3] This dependence should theoretically 
allow one to place limits upon the WIMP mass in the event of ei- 
ther a positive or null detection of dark stars at the Galactic centre. 



6.3 Binaries and higher-multiplicity stellar systems 

Although we have not calculated the WIMP luminosities achiev- 
able by systems consisting of more than one star, our results with 
single stars make it easy to comment on this scenario. If a bi- 
nary system's internal orbital plane was partially aligned with its 
gross orbit about the Galactic centre, the motions of its component 
stars would at various stages counteract the motion of the system 
about the central black hole. At certain times the component stars 
would have far smaller velocities relative to the WIMP halo than 
if they were orbiting as single stars on a similar orbit, resulting in 
increased capture rates. The systems most effective at achieving 
this boost would be those with the highest orbital speeds, as these 
would produce the greatest reduction in the relative velocity be- 
tween stars and WIMPs. This effect would therefore be strongest 
in short-period, low-separation, higher-mass binary systems. Since 
the effects of WIMP annihilation are most marked in low-mass 
stars though, the optimal configuration would be a close binary 
with a maximal mass difference. A system consisting of ~1Mq 
and ~4Mq partners orbiting with a period of 5 hr would have an 
orbital velocity of ~700 km s~ 1 , enough to have a profound impact 
upon capture rates shown in Fig. [TT] for example. In favourable 
cases, binary systems could mimic the effects of highly eccentric 
orbits upon capture rates, allowing stars on almost circular orbits to 
achieve significant capture rates, and further boosting capture rates 
from elliptical orbits. Similar effects could be expected in systems 
consisting of three or more stars, though their stability on orbits 
close to the Galactic Centre might be doubtful; indeed, even bina- 
ries might not survive for long near the central black hole l |Perets| 
|2009l >. 

6.4 Observational constraints and prospects 

Single stars on circular orbits capturing dark matter from an isother- 
mal halo cannot achieve WIMP luminosities any greater than 1 
per cent of their nuclear luminosities. Even if the WIMP velocities 
were not isothermal, but instead followed the A-body distribution, 
Fig. [T3] indicates that capture would not be boosted by more than 
an additional factor of 5. From the results of Sect.[3] we know that 
log 10 [Lw,max/inuc(0)] < — 1 would not result in any significant 
change to a star's structure or evolution. Whilst this level of WIMP 
luminosity would create small convective cores in some stars, po- 
tentially interesting for asteroseismology of Galactic centre popu- 
lations some time in the distant future, we find it very unlikely that 
main-sequence dark stars would exist outside binaries on any cir- 
cular orbits in the Milky Way. Stars following elliptical orbits are 
not only far more likely to be dark stars, but are also considerably 
more common at the Galactic centre. 

The central 30 pc of the Milky Way not only contains a 
3 — 4 x 10 6 Mq black hole, but also two of the densest star clusters 
in the Galaxy, including the Arches cluster. In the late 1980s, an un- 
usual star with broad H I and He I emission lines was detected less 
than 0.5 pc from the central compact radio source, SgrA*, which 
is thought to be associated with the central black hole ( [Forrest| 
|et al.||1986| |Allen et al.|1990) . Because the centre of the Galaxy 
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is shrouded with dust, observations can only take place in the in- 
frared, so the normal spectral information used to identify stars is 
not available. Over the next few years, an increasing number of 
such stars were discovered, appearing to be helium-rich blue super- 
giants and Wolf-Rayet stars, with masses of up to 100 Mo (Krabbe 
|et al.|199"T"[ |1995| >. The presence of such young stars close to the 
Galactic centre was difficult to understand, as it was unclear how a 
gas cloud could condense to form a star in a region so close to the 
central black hole due to the extreme tidal forces there. This prob- 
lem came to be known as the 'Paradox of Youth' (Sander s)l992| 
|Morris|1993| >. Recent work has increased the number of known OB 
stars in the central parsec of the Galaxy to close to 100, excluding 
the central square arcsecond (Paumard et al. 2006). These young 
stars appear to form two counter rotating discs, suggesting that they 
are possibly associated with different star formation events in dense 
accreting matter (|Levin & B eloborodo vp003| |Genzel et al.|2003| 
|Paumard et al.|2006) . 

The very fact that we know the mass of the central black hole 
itself is due to the discovery and subsequent tracking of stars even 
closer than the young He I stars: a cluster of stars was discovered 
in the 1990s within the central square arcsecond, and dubbed the 
'SgrA* stellar cluster' (Ge nzel et al.|1997) . The stars in this cluster 
(referred to as S stars, E stars or SO stars apparently depending 
upon the native language of the lead researcher) move extremely 
close to the central black hole. The star which has been observed 
moving closest to the central black hole is called S14, E2 or SO-16, 
and was seen within 45 AU, or only 600 Schwarszchild radii, of the 
black hole (Ghez et al. 2005 ). The kinematics of the SgrA* cluster 
is such that stars are on randomly oriented, highly elliptical orbits, 
rather than circular ones confined to a single disc. 

Near-infrared observations of the S stars made with the Keck 
telescope and VLT revealed that their atmospheres did not contain 
CO ( |Genzel et al.||199"7} , setting a lower bound on their surface 
temperatures. Further work led to the conclusion that the S stars are 
in fact simply 10-15 Mq main-sequence stars (Ghez et al .|2005| 
Martins et al. 2008). The short main-sequence lifetimes of such 
high-mass stars (~10 7 yr) presents a further Paradox of Youth at the 
Galactic centre, not explainable by star formation in an accretion 
disc due to the random orientation of the S-star orbits. 

Many have tried to explain the presence of the S stars. One 
idea is that they formed far from the Galactic centre (perhaps in 
the Arches cluster) and subsequently migrated inwards ( |Gerhard| 
|2001[ >. Another is that they formed in situ during an earlier era when 
the density was much larger than it is today (Levin & Belo borodov| 
|2003| ). They could also be old stars which look young because they 
have collided with other stars (Genzel et al. 2003), a scenario rem- 
iniscent of the blue straggler phenomenon in globular clusters. A 
rather convincing explanation is that they were originally members 
of binaries belonging to one of the outer discs, which were per- 
turbed either by interactions between the discs (Lockmann et al 



2008 



or by interactions with other massive objects (Perets et al. 



2007[ l. The three-body interaction then caused one star from each 
binary to become tightly bound to the black hole, and the other to 
be ejected as a hypervelocity star. 

What implications might dark stars have for this picture? We 
have already seen that stars burning dark matter have significantly 
increased main-sequence lifetimes. One might then imagine a sce- 
nario whereby stars are created elsewhere and migrate to the Galac- 
tic centre, where the presence of dark matter extends their lifetimes. 
This might provide an alternative explanation for either the S stars 
or the outer stellar discs of OB-type stars. However, such an ex- 
planation is incomplete. The problem with models where stars are 



created elsewhere and migrate to the centre of the Galaxy is that the 
inspiralling timescale is typically very large compared to their main 
sequence lifetime. One would therefore expect that stars should 
have left the main sequence by the time they arrive at the central 
region. Furthermore, we have shown that more massive stars re- 
quire higher dark matter densities than low mass ones to experience 
any structural changes; it is highly unlikely that a star as massive as 
10 M© could capture enough WIMPs to significantly alter its main- 
sequence evolution on any realistic orbit near the Galactic centre. 
One possibility is that such a star could reach the end of its main 
sequence lifetime during the migration, arrive at the Galactic cen- 
tre and then begin capturing large numbers of WIMPs. If burning 
WIMPs during its post-main-sequence evolution made such a star 
begin to resemble an OB or Wolf-Rayet star, or revert to looking 
outwardly like a main sequence star, this could provide an addi- 
tional explanation for the dense stellar discs or the S stars, respec- 
tively. We will consider the prospect of post-main-sequence dark 
stars in a later paper. 

Whilst such an explanation for the Paradox of Youth must be 
considered improbable, it is interesting to remember that the S stars 
are indeed on more elliptical orbits than other stars at the Galactic 
centre, which would be consistent with them having accreted far 
more dark matter than others (Zh u et al.|2008| l. 

More promising is the possibility that future observations of 
the Galactic centre will reveal fainter, lower-mass stars (although 
there is some suggestion that the initial mass function within the 
central parsec could be top heavy; Maness et al. 2007). If the bi- 
nary disruption scenario is indeed the source of the S stars, one 
would expect that the bursts of star formation which created them 
in the outer discs would also have produced lower-mass stars. Some 
such stars would form in binaries, and could conceivably follow 
the same path to the Galactic centre as the S stars. This would pro- 
duce a population of low-mass stars in the central square arcsecond 
with randomly-oriented, highly elliptical orbits similar to those of 
the S stars. Likewise, most of the other explanations for the origin 
of the S stars could also involve formation and subsequent migra- 
tion/disruption of a lower-mass population alongside the S stars, 
also leading to a population of potential low-mass dark stars. 

Finally, the prospect of dark stars forming in binary systems 
opens a very promising channel through which they might be ob- 
served. Stars within a binary can be compared photometrically if 
they have similar brightness, allowing their masses and evolution- 
ary states to be determined. In some cases, binaries might consist 
of a low-mass star which is significantly affected by WIMP cap- 
ture and annihilation, and a high-mass partner which is too mas- 
sive to show any effects whatsoever. The most striking example of 
this would be a binary consisting of a low-mass star 'frozen' by 
WIMP burning (resembling a protostar) and a higher-mass com- 
panion which had evolved all the way into a white dwarf, though 
such a system would be difficult to observe at the Galactic centre 
because of the faintness of the white dwarf. 



7 CONCLUSIONS 

When the energy injected due to the annihilation of WIMPs ap- 
proaches that of nuclear burning, the capture of weakly interacting 
dark matter will significantly alter the structure and evolution of 
stars on the main sequence. Stars on circular orbits in the Milky 
Way are extremely unlikely to achieve sufficient capture rates for 
this to occur unless they are present in binaries. Stars orbiting close 
to the Galactic centre on elliptical orbits have their capture rates 
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strongly boosted in comparison to those on circular orbits. The 
velocity distribution of dark matter near the Galactic centre may 
be highly non-Gaussian, further boosting capture rates on ellip- 
tical orbits by nearly an order of magnitude. Assuming that the 
nuclear-scattering cross-sections are equal to their current exper- 
imental limits, that dark matter forms a spike around the super- 
massive black hole at the Galactic centre, and that the dark matter 
distribution on larger scales has undergone adiabatic contraction, 
stars of 1 M Q and below will reach break-even between annihila- 
tion and fusion energy on orbits with periods of up to 50 years and 
eccentricities as low as 0.9. 1.5 M Q stars can achieve the same goal 
with comparable orbital eccentricities if they orbit the central black 
hole in 10 years or less. Without adiabatic contraction of the galac- 
tic halo, orbits at least as short as this and eccentricities of about 
0.99 are required for stars of a solar mass and below to become 
dark stars. 

These requirements are likely to be significantly relaxed for 
stars in binary systems. A binary consisting of a low-mass protostar 
and a highly-evolved massive star would make the impact of WIMP 
annihilation very hard to deny. 

The observation of one or more stars at the Galactic centre ex- 
hibiting the properties we have described would strongly suggest 
the influence of WIMP dark matter. Conversely, since we have as- 
sumed scattering cross-sections compatible with the current exper- 
imental limits, the observation of even a single completely normal 
star or binary on the orbits we have discussed would allow one to 
place stringent limits on the properties of dark matter and its den- 
sity at the Galactic centre. If one instead assumed a particular halo 
model for dark matter at the Galactic centre, the dependence of the 
capture rate upon the WIMP mass and spin-dependent scattering 
cross-section would allow one to derive limits on these parameters 
which are highly competitive with current direct-detection sensitiv- 
ities. If a star were seen on an orbit where we expect effects even 
without adiabatic contraction of the Galactic halo on large scales 
(i.e. M < 1M Q , P < lOyr, e > 0.99), then the derived limits 
could be made mostly independent of the halo model. 
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